library(MASS)
library(margins)
library(clubSandwich)
library(stargazer)
library(survival)



data<-read.csv("datafull_190123.csv")
data<-data[-c(1)]

### Summary Statistics (Appendix, Table 1)

sumstats<-data[,c("polity2","regimechange","democracy","death","agepre","agepost","tenpre","tenpost","gdppcln","personal","otherirreg","allcoup_bin")]
#stargazer(sumstats,out="sumstats", label="tab:sumstats",title="Summary Statistics")


### H0 Analysis - effect of leader death on different measures of political change.  

#DV: Polity

#1-year DV (Appendix, Table 2)

nondem<-subset(data,polity_lag1<0)
nondem<-subset(nondem, is.na(death)==FALSE )

d1<-subset(nondem,death_for1==0 )
d1<-subset(d1,is.na(polity_for1)==FALSE)

d<-d1
m1<- lm(polity_for1~death,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE)
m2<- lm(polity_for1~death+polity_lag1,data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m2)
pol1_noage_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m3<- lm(polity_for1~death+agepre+polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m3)
pol1_age_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(polity_for1) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m4<- lm(polity_for1~death+agepre+polity_lag1,data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m5<- lm(polity_for1~death+agepre+tenpre+polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))


d1<-subset(d1,is.na(ccode)==FALSE &is.na(year)==FALSE)
d<-d1
m6<- lm(polity_for1~death+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE)
m7<- lm(polity_for1~death+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m7,cluster=d$ccode,type="CR1")
se7<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m7)
pol1_noage_fe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m8<- lm(polity_for1~death+agepre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m8,cluster=d$ccode,type="CR1")
se8<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m8)
pol1_age_fe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(polity_for1) ==FALSE & is.na(ccode)==FALSE &is.na(year)==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m9<- lm(polity_for1~death+agepre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m9,cluster=d$ccode,type="CR1")
se9<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m10<- lm(polity_for1~death+agepre+tenpre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m10,cluster=d$ccode,type="CR1")
se10<-sqrt(diag(vc))

plotty_pol1<-rbind(pol1_noage_nofe,pol1_age_nofe,pol1_noage_fe,pol1_age_fe)



#5-year DV (Appendix, Table 3)


nondem<-subset(data,polity_lag1<0)
nondem<-subset(nondem, is.na(death)==FALSE )

d5<-subset(nondem,death_for5==0 )
d5<-subset(d5,is.na(polity_for5)==FALSE)

d<-d5
m1<- lm(polity_for5~death,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d<-subset(d5,is.na(polity_lag1)==FALSE)
m2<- lm(polity_for5~death+polity_lag1,data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m2)
pol5_noage_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d5,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m3<- lm(polity_for5~death+agepre+polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m3)
pol5_age_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(polity_for5) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m4<- lm(polity_for5~death+agepre+polity_lag1,data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-subset(d5,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m5<- lm(polity_for5~death+agepre+tenpre+polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))


d5<-subset(d5,is.na(ccode)==FALSE &is.na(year)==FALSE)

d<-d5
m6<- lm(polity_for5~death+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))

d<-subset(d5,is.na(polity_lag1)==FALSE)
m7<- lm(polity_for5~death+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m7,cluster=d$ccode,type="CR1")
se7<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m7)
pol5_noage_fe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d5,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m8<- lm(polity_for5~death+agepre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m8,cluster=d$ccode,type="CR1")
se8<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m8)
pol5_age_fe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(polity_for5) ==FALSE & is.na(ccode)==FALSE &is.na(year)==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m9<- lm(polity_for5~death+agepre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m9,cluster=d$ccode,type="CR1")
se9<-sqrt(diag(vc))

d<-subset(d5,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m10<- lm(polity_for5~death+agepre+tenpre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m10,cluster=d$ccode,type="CR1")
se10<-sqrt(diag(vc))

plotty_pol5<-rbind(pol5_noage_nofe,pol5_age_nofe,pol5_noage_fe,pol5_age_fe)



#10-year DV (Appendix, Table 4)


nondem<-subset(data,polity_lag1<0)
nondem<-subset(nondem, is.na(death)==FALSE )

d10<-subset(nondem,death_for10==0 )
d10<-subset(d10,is.na(polity_for10)==FALSE)

d<-d10
m1<- lm(polity_for10~death,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d<-subset(d10,is.na(polity_lag1)==FALSE)
m2<- lm(polity_for10~death+polity_lag1,data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m2)
pol10_noage_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d10,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m3<- lm(polity_for10~death+agepre+polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m3)
pol10_age_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(polity_for10) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m4<- lm(polity_for10~death+agepre+polity_lag1,data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-subset(d10,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m5<- lm(polity_for10~death+agepre+tenpre+polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))


d10<-subset(d10,is.na(ccode)==FALSE &is.na(year)==FALSE)

d<-d10
m6<- lm(polity_for10~death+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))

d<-subset(d10,is.na(polity_lag1)==FALSE)
m7<- lm(polity_for10~death+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m7,cluster=d$ccode,type="CR1")
se7<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m7)
pol10_noage_fe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d10,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m8<- lm(polity_for10~death+agepre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m8,cluster=d$ccode,type="CR1")
se8<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m8)
pol10_age_fe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(polity_for10) ==FALSE & is.na(ccode)==FALSE &is.na(year)==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m9<- lm(polity_for10~death+agepre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m9,cluster=d$ccode,type="CR1")
se9<-sqrt(diag(vc))

d<-subset(d10,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m10<- lm(polity_for10~death+agepre+tenpre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m10,cluster=d$ccode,type="CR1")
se10<-sqrt(diag(vc))

plotty_pol10<-rbind(pol10_noage_nofe,pol10_age_nofe,pol10_noage_fe,pol10_age_fe)




## DV: Democratization 


#1-year DV (Appendix, Table 5)

nondem<-subset(data,democracy_lag1==0)
nondem<-subset(nondem, is.na(death)==FALSE )

d1<-subset(nondem,death_for1==0 )
d1<-subset(d1,is.na(democracy_for1)==FALSE)

d<-d1
m1<- glm(democracy_for1~death,data=d,family=binomial(link=logit))
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE)
m2<- glm(democracy_for1~death+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m2)
dem1_noage_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m3<- glm(democracy_for1~death+agepre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m3)
dem1_age_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(democracy_for1) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m4<- glm(democracy_for1~death+agepre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m5<- glm(democracy_for1~death+agepre+tenpre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))


d1<-subset(d1,is.na(ccode)==FALSE &is.na(year)==FALSE)

d<-d1
m6<- clogit(democracy_for1~death+strata(ccode),data=d)

d<-subset(d1,is.na(polity_lag1)==FALSE)
m7<- clogit(democracy_for1~death+polity_lag1+strata(ccode),data=d)
coef<-coefficients(m7)[1]
se<-coef(summary(m7))[1,3]
dem1_noage_fe<-c(coef,coef+1.96*se,coef-1.96*se)

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m8<- clogit(democracy_for1~death+agepre+polity_lag1+strata(ccode),data=d)
coef<-coefficients(m8)[1]
se<-coef(summary(m8))[1,3]
dem1_age_fe<-c(coef,coef+1.96*se,coef-1.96*se)

d<-subset(nondem,is.na(democracy_for1) ==FALSE & is.na(ccode)==FALSE &is.na(year)==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m9<- clogit(democracy_for1~death+agepre+polity_lag1+strata(ccode),data=d)


d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m10<- clogit(democracy_for1~death+agepre+tenpre+polity_lag1+strata(ccode),data=d)


plotty_dem1<-rbind(dem1_noage_nofe,dem1_age_nofe,dem1_noage_fe,dem1_age_fe)


## 5-year DV (Appendix, Table 6)

nondem<-subset(data,democracy_lag1==0)
nondem<-subset(nondem, is.na(death)==FALSE )

d5<-subset(nondem,death_for5==0 )
d5<-subset(d5,is.na(democracy_for5)==FALSE)

d<-d5
m1<- glm(democracy_for5~death,data=d,family=binomial(link=logit))
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d<-subset(d5,is.na(polity_lag1)==FALSE)
m2<- glm(democracy_for5~death+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m2)
dem5_noage_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d5,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m3<- glm(democracy_for5~death+agepre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m3)
dem5_age_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(democracy_for5) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m4<- glm(democracy_for5~death+agepre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-subset(d5,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m5<- glm(democracy_for5~death+agepre+tenpre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))


d5<-subset(d5,is.na(ccode)==FALSE &is.na(ccode)==FALSE)

d<-d5
m6<- clogit(democracy_for5~death+strata(ccode),data=d)

d<-subset(d5,is.na(polity_lag1)==FALSE)
m7<- clogit(democracy_for5~death+polity_lag1+strata(ccode),data=d)
coef<-coefficients(m7)[1]
se<-coef(summary(m7))[1,3]
dem5_noage_fe<-c(coef,coef+1.96*se,coef-1.96*se)

d<-subset(d5,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m8<- clogit(democracy_for5~death+agepre+polity_lag1+strata(ccode),data=d)
coef<-coefficients(m8)[1]
se<-coef(summary(m8))[1,3]
dem5_age_fe<-c(coef,coef+1.96*se,coef-1.96*se)


d<-subset(nondem,is.na(democracy_for5) ==FALSE & is.na(ccode)==FALSE &is.na(year)==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m9<- clogit(democracy_for5~death+agepre+polity_lag1+strata(ccode),data=d)

d<-subset(d5,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m10<- clogit(democracy_for5~death+agepre+tenpre+polity_lag1+strata(ccode),data=d)

plotty_dem5<-rbind(dem5_noage_nofe,dem5_age_nofe,dem5_noage_fe,dem5_age_fe)


## 10-year DV (Appendix, Table 7)

nondem<-subset(data,democracy_lag1==0)
nondem<-subset(nondem, is.na(death)==FALSE )

d10<-subset(nondem,death_for10==0 )
d10<-subset(d10,is.na(democracy_for10)==FALSE)

d<-d10
m1<- glm(democracy_for10~death,data=d,family=binomial(link=logit))
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d<-subset(d10,is.na(polity_lag1)==FALSE)
m2<- glm(democracy_for10~death+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m2)
dem10_noage_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d10,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m3<- glm(democracy_for10~death+agepre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m3)
dem10_age_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(democracy_for10) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m4<- glm(democracy_for10~death+agepre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-subset(d10,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m5<- glm(democracy_for10~death+agepre+tenpre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))


d5<-subset(d10,is.na(ccode)==FALSE &is.na(ccode)==FALSE)

d<-d10
m6<- clogit(democracy_for10~death+strata(ccode),data=d)

d<-subset(d10,is.na(polity_lag1)==FALSE)
m7<- clogit(democracy_for10~death+polity_lag1+strata(ccode),data=d)
coef<-coefficients(m7)[1]
se<-coef(summary(m7))[1,3]
dem10_noage_fe<-c(coef,coef+1.96*se,coef-1.96*se)

d<-subset(d10,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m8<- clogit(democracy_for10~death+agepre+polity_lag1+strata(ccode),data=d)
coef<-coefficients(m8)[1]
se<-coef(summary(m8))[1,3]
dem10_age_fe<-c(coef,coef+1.96*se,coef-1.96*se)


d<-subset(nondem,is.na(democracy_for10) ==FALSE & is.na(ccode)==FALSE &is.na(year)==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m9<- clogit(democracy_for10~death+agepre+polity_lag1+strata(ccode),data=d)

d<-subset(d10,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m10<- clogit(democracy_for10~death+agepre+tenpre+polity_lag1+strata(ccode),data=d)

plotty_dem10<-rbind(dem10_noage_nofe,dem10_age_nofe,dem10_noage_fe,dem10_age_fe)


### Regime Change

## 1-Year DV (Appendix, Table 8)

nondem<-subset(data,military_lag1==1 |monarchy_lag1==1|personal_lag1==1|party_lag1==1)
nondem<-subset(nondem, is.na(death)==FALSE )

d1<-subset(nondem,death_for1==0 )
d1<-subset(d1,is.na(regimechange_for1)==FALSE)

d<-d1
m1<- glm(regimechange_for1~death,data=d,family=binomial(link=logit))
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE)
m2<- glm(regimechange_for1~death+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m2)
rc1_noage_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m3<- glm(regimechange_for1~death+agepre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m3)
rc1_age_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(regimechange_for1) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m4<- glm(regimechange_for1~death+agepre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m5<- glm(regimechange_for1~death+agepre+tenpre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))


d1<-subset(d1,is.na(ccode)==FALSE &is.na(year)==FALSE)

d<-d1
m6<- clogit(regimechange_for1~death+strata(ccode),data=d)

d<-subset(d1,is.na(polity_lag1)==FALSE)
m7<- clogit(regimechange_for1~death+polity_lag1+strata(ccode),data=d)
coef<-coefficients(m7)[1]
se<-coef(summary(m7))[1,3]
rc1_noage_fe<-c(coef,coef+1.96*se,coef-1.96*se)

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m8<- clogit(regimechange_for1~death+agepre+polity_lag1+strata(ccode),data=d)
coef<-coefficients(m8)[1]
se<-coef(summary(m8))[1,3]
rc1_age_fe<-c(coef,coef+1.96*se,coef-1.96*se)

d<-subset(nondem,is.na(regimechange_for1) ==FALSE & is.na(ccode)==FALSE &is.na(year)==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m9<- clogit(democracy_for1~death+agepre+polity_lag1+strata(ccode),data=d)


d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m10<- clogit(regimechange_for1~death+agepre+tenpre+polity_lag1+strata(ccode),data=d)


plotty_rc1<-rbind(rc1_noage_nofe,rc1_age_nofe,rc1_noage_fe,rc1_age_fe)


## 5-Yr DV (Appendix, Table 5)

nondem<-subset(data,military_lag1==1 |monarchy_lag1==1|personal_lag1==1|party_lag1==1)
nondem<-subset(nondem, is.na(death)==FALSE )

d5<-subset(nondem,death_for5==0 )
d5<-subset(d5,is.na(regimechange_for5)==FALSE)

d<-d5
m1<- glm(regimechange_for5~death,data=d,family=binomial(link=logit))
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d<-subset(d5,is.na(polity_lag1)==FALSE)
m2<- glm(regimechange_for5~death+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m2)
rc5_noage_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d5,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m3<- glm(regimechange_for5~death+agepre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m3)
rc5_age_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(regimechange_for5) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m4<- glm(regimechange_for5~death+agepre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m5<- glm(regimechange_for5~death+agepre+tenpre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))

d5<-subset(d5,is.na(ccode)==FALSE &is.na(year)==FALSE)
d<-d5

m6<- clogit(regimechange_for5~death+strata(ccode),data=d)

d<-subset(d5,is.na(polity_lag1)==FALSE)
m7<- clogit(regimechange_for5~death+polity_lag1+strata(ccode),data=d)
coef<-coefficients(m7)[1]
se<-coef(summary(m7))[1,3]
rc5_noage_fe<-c(coef,coef+1.96*se,coef-1.96*se)

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m8<- clogit(regimechange_for5~death+agepre+polity_lag1+strata(ccode),data=d)
coef<-coefficients(m8)[1]
se<-coef(summary(m8))[1,3]
rc5_age_fe<-c(coef,coef+1.96*se,coef-1.96*se)

d<-subset(nondem,is.na(regimechange_for5) ==FALSE & is.na(ccode)==FALSE &is.na(year)==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m9<- clogit(democracy_for5~death+agepre+polity_lag1+strata(ccode),data=d)


d<-subset(d5,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m10<- clogit(regimechange_for5~death+agepre+tenpre+polity_lag1+strata(ccode),data=d)


plotty_rc5<-rbind(rc5_noage_nofe,rc5_age_nofe,rc5_noage_fe,rc5_age_fe)


## 10-yr DV (Appendix, Table 10)


nondem<-subset(data,military_lag1==1 |monarchy_lag1==1|personal_lag1==1|party_lag1==1)
nondem<-subset(nondem, is.na(death)==FALSE )

d10<-subset(nondem,death_for10==0 )
d10<-subset(d10,is.na(regimechange_for10)==FALSE)

d<-d10
m1<- glm(regimechange_for10~death,data=d,family=binomial(link=logit))
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d<-subset(d10,is.na(polity_lag1)==FALSE)
m2<- glm(regimechange_for10~death+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m2)
rc10_noage_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d10,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m3<- glm(regimechange_for10~death+agepre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m3)
rc10_age_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(regimechange_for10) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m4<- glm(regimechange_for10~death+agepre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m5<- glm(regimechange_for10~death+agepre+tenpre+polity_lag1,data=d,family=binomial(link=logit))
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))

d5<-subset(d10,is.na(ccode)==FALSE &is.na(year)==FALSE)
d<-d10

m6<- clogit(regimechange_for10~death+strata(ccode),data=d)

d<-subset(d10,is.na(polity_lag1)==FALSE)
m7<- clogit(regimechange_for10~death+polity_lag1+strata(ccode),data=d)
coef<-coefficients(m7)[1]
se<-coef(summary(m7))[1,3]
rc10_noage_fe<-c(coef,coef+1.96*se,coef-1.96*se)

d<-subset(d10,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m8<- clogit(regimechange_for10~death+agepre+polity_lag1+strata(ccode),data=d)
coef<-coefficients(m8)[1]
se<-coef(summary(m8))[1,3]
rc10_age_fe<-c(coef,coef+1.96*se,coef-1.96*se)

d<-subset(nondem,is.na(regimechange_for10) ==FALSE & is.na(ccode)==FALSE &is.na(year)==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m9<- clogit(democracy_for10~death+agepre+polity_lag1+strata(ccode),data=d)


d<-subset(d10,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m10<- clogit(regimechange_for10~death+agepre+tenpre+polity_lag1+strata(ccode),data=d)


plotty_rc10<-rbind(rc10_noage_nofe,rc10_age_nofe,rc10_noage_fe,rc10_age_fe)


### Absolute change in polity


#1-year DV (Appendix, Table 11)

nondem<-subset(data,polity_lag1<0)
nondem<-subset(nondem, is.na(death)==FALSE )

d1<-subset(nondem,death_for1==0 )
d1<-subset(d1,is.na(abspolchange_for1)==FALSE)

d<-d1
m1<- lm(abspolchange_for1~death,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE)
m2<- lm(abspolchange_for1~death+polity_lag1,data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m2)
abs1_noage_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m3<- lm(abspolchange_for1~death+agepre+polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m3)
abs1_age_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(abspolchange_for1) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m4<- lm(abspolchange_for1~death+agepre+polity_lag1,data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m5<- lm(abspolchange_for1~death+agepre+tenpre+polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))

d1<-subset(d1,is.na(ccode)==FALSE &is.na(year)==FALSE)
d<-d1
m6<- lm(abspolchange_for1~death+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE)
m7<- lm(abspolchange_for1~death+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m7,cluster=d$ccode,type="CR1")
se7<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m7)
abs1_noage_fe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m8<- lm(abspolchange_for1~death+agepre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m8,cluster=d$ccode,type="CR1")
se8<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m8)
abs1_age_fe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(polity_for1) ==FALSE & is.na(ccode)==FALSE &is.na(year)==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m9<- lm(abspolchange_for1~death+agepre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m9,cluster=d$ccode,type="CR1")
se9<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m10<- lm(abspolchange_for1~death+agepre+tenpre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m10,cluster=d$ccode,type="CR1")
se10<-sqrt(diag(vc))

plotty_abs1<-rbind(abs1_noage_nofe,abs1_age_nofe,abs1_noage_fe,abs1_age_fe)



#5-year DV (Appendix, Table 12)

nondem<-subset(data,polity_lag1<0)
nondem<-subset(nondem, is.na(death)==FALSE )

d5<-subset(nondem,death_for5==0 )
d5<-subset(d5,is.na(abspolchange_for5)==FALSE)

d<-d5
m1<- lm(abspolchange_for5~death,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d<-subset(d5,is.na(polity_lag1)==FALSE)
m2<- lm(abspolchange_for5~death+polity_lag1,data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m2)
abs5_noage_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d5,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m3<- lm(abspolchange_for5~death+agepre+polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m3)
abs5_age_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(abspolchange_for5) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m4<- lm(abspolchange_for5~death+agepre+polity_lag1,data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-subset(d5,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m5<- lm(abspolchange_for5~death+agepre+tenpre+polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))

d5<-subset(d1,is.na(ccode)==FALSE &is.na(year)==FALSE)
d<-d5
m6<- lm(abspolchange_for5~death+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))

d<-subset(d5,is.na(polity_lag1)==FALSE)
m7<- lm(abspolchange_for5~death+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m7,cluster=d$ccode,type="CR1")
se7<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m7)
abs5_noage_fe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d5,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m8<- lm(abspolchange_for5~death+agepre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m8,cluster=d$ccode,type="CR1")
se8<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m8)
abs5_age_fe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(polity_for5) ==FALSE & is.na(ccode)==FALSE &is.na(year)==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m9<- lm(abspolchange_for5~death+agepre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m9,cluster=d$ccode,type="CR1")
se9<-sqrt(diag(vc))

d<-subset(d5,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m10<- lm(abspolchange_for5~death+agepre+tenpre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m10,cluster=d$ccode,type="CR1")
se10<-sqrt(diag(vc))

plotty_abs5<-rbind(abs5_noage_nofe,abs5_age_nofe,abs5_noage_fe,abs5_age_fe)


## 10 year DV (Appendix, Table 13)

nondem<-subset(data,polity_lag1<0)
nondem<-subset(nondem, is.na(death)==FALSE )

d10<-subset(nondem,death_for5==0 )
d10<-subset(d10,is.na(abspolchange_for10)==FALSE)

d<-d10
m1<- lm(abspolchange_for10~death,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d<-subset(d10,is.na(polity_lag1)==FALSE)
m2<- lm(abspolchange_for10~death+polity_lag1,data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m2)
abs10_noage_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d10,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m3<- lm(abspolchange_for10~death+agepre+polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m3)
abs10_age_nofe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(abspolchange_for10) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m4<- lm(abspolchange_for10~death+agepre+polity_lag1,data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-subset(d10,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m5<- lm(abspolchange_for10~death+agepre+tenpre+polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))

d10<-subset(d10,is.na(ccode)==FALSE &is.na(year)==FALSE)
d<-d10
m6<- lm(abspolchange_for10~death+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))

d<-subset(d10,is.na(polity_lag1)==FALSE)
m7<- lm(abspolchange_for10~death+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m7,cluster=d$ccode,type="CR1")
se7<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m7)
abs10_noage_fe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(d10,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m8<- lm(abspolchange_for10~death+agepre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m8,cluster=d$ccode,type="CR1")
se8<-sqrt(diag(vc))
se<-sqrt(diag(vc))
coef<-coefficients(m8)
abs10_age_fe<-c(coef[2],coef[2]+1.96*se[2],coef[2]-1.96*se[2])

d<-subset(nondem,is.na(polity_for10) ==FALSE & is.na(ccode)==FALSE &is.na(year)==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m9<- lm(abspolchange_for10~death+agepre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m9,cluster=d$ccode,type="CR1")
se9<-sqrt(diag(vc))

d<-subset(d10,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m10<- lm(abspolchange_for10~death+agepre+tenpre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m10,cluster=d$ccode,type="CR1")
se10<-sqrt(diag(vc))

plotty_abs10<-rbind(abs10_noage_nofe,abs10_age_nofe,abs10_noage_fe,abs10_age_fe)


###H0 Plot - Figure 2 in the main article. 

pdf(file="figH0.pdf")
par(mfrow=c(1,1))
colors<-c("gray75", "gray50","gray25","gray0")

plot_1<-plotty_pol1
plot_5<-plotty_pol5
plot_10<-plotty_pol10


plot(x=plot_1[,1],y=c(73:70),xlim=c(-5,5),ylim=c(1,73),pch=20,col=colors,yaxt="n",xlab="Coefficient Estimate",ylab="",cex=0.8)
segments(x0=plot_1[,2],y0=c(73:70),x1=plot_1[,3],y1=c(73:70),col=colors)
points(x=plot_5[,1],y=c(67:64),pch=20,col=colors)
segments(x0=plot_5[,2],y0=c(67:64),x1=plot_5[,3],y1=c(67:64),col=colors)
points(x=plot_10[,1],y=c(61:58),pch=20,col=colors)
segments(x0=plot_10[,2],y0=c(61:58),x1=plot_10[,3],y1=c(61:58),col=colors)
abline(h=56)

plot_1<-plotty_dem1
plot_5<-plotty_dem5
plot_10<-plotty_dem10
points(x=plot_1[,1],y=c(54:51),pch=20,col=colors)
segments(x0=plot_1[,2],y0=c(54:51),x1=plot_1[,3],y1=c(54:51),col=colors)
points(x=plot_5[,1],y=c(48:45),pch=20,col=colors)
segments(x0=plot_5[,2],y0=c(48:45),x1=plot_5[,3],y1=c(48:45),col=colors)
points(x=plot_10[,1],y=c(42:39),pch=20,col=colors)
segments(x0=plot_10[,2],y0=c(42:39),x1=plot_10[,3],y1=c(42:39),col=colors)
abline(h=37)

plot_1<-plotty_rc1
plot_5<-plotty_rc5
plot_10<-plotty_rc10
colors<-c("gray75", "gray50","gray25","gray0")
points(x=plot_1[,1],y=c(35:32),pch=20,col=colors)
segments(x0=plot_1[,2],y0=c(35:32),x1=plot_1[,3],y1=c(35:32),col=colors)
points(x=plot_5[,1],y=c(29:26),pch=20,col=colors)
segments(x0=plot_5[,2],y0=c(29:26),x1=plot_5[,3],y1=c(29:26),col=colors)
points(x=plot_10[,1],y=c(23:20),pch=20,col=colors)
segments(x0=plot_10[,2],y0=c(23:20),x1=plot_10[,3],y1=c(23:20),col=colors)
abline(h=18)

plot_1<-plotty_abs1
plot_5<-plotty_abs5
plot_10<-plotty_abs10
points(x=plot_1[,1],y=c(16:13),pch=20,col=colors)
segments(x0=plot_1[,2],y0=c(16:13),x1=plot_1[,3],y1=c(16:13),col=colors)
points(x=plot_5[,1],y=c(10:7),pch=20,col=colors)
segments(x0=plot_5[,2],y0=c(10:7),x1=plot_5[,3],y1=c(10:7),col=colors)
points(x=plot_10[,1],y=c(4:1),pch=20,col=colors)
segments(x0=plot_10[,2],y0=c(4:1),x1=plot_10[,3],y1=c(4:1),col=colors)

abline(v=0,lty=2)

axis(2, at=c(2.5,7.5,12.5,21.5,27.5,33.5,40.5,46.5,52.5,59.5,65.5,71.5),tick=FALSE,labels=rep(c("10yr","5yr","1yr"),4), las=2,cex=0.8)
legend(2.75,72.5,c("No age, no FE","Age, no FE","No age, FE","Age, FE"), pch=20,col=colors,text.font=3,cex=0.8)
abline(v=0,lty=2)

text(-5.25,59.5,"DV: Polity",vfont=c("sans serif","bold italic"),pos=4)

text(-5.25,40.5,"DV: Democracy",vfont=c("sans serif","bold italic"),pos=4)
text(-5.25,21.5,"DV: Regime Change",vfont=c("sans serif","bold italic"),pos=4)
text(-5.25,3.5,"DV: Absolute",vfont=c("sans serif","bold italic"),pos=4)
text(-5.15,0,"Polity Change",vfont=c("sans serif","bold italic"),pos=4)
dev.off()


####### H1 Analysis:  Effect of leader death and age on political liberalization ########

# 1-Year DV (Appendix, Table 14)

nondem<-subset(data,polity_lag1<0)
nondem<-subset(nondem, is.na(death)==FALSE & is.na(agepre)==FALSE )

d1<-subset(nondem,death_for1==0 )
d1<-subset(d1,is.na(polity_for1)==FALSE)

d<-d1
m1<- lm(polity_for1~death*agepre,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE)
m2<- lm(polity_for1~death*agepre + polity_lag1,data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))
a<-seq(quantile(d$agepre,0.1, na.rm=TRUE),quantile(d$agepre,0.9,na.rm=TRUE),1)
marg1_nofe<-margins(m2,at=list(agepre=a),vcov=vc,variables="death")

d<-subset(nondem,is.na(polity_for1) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m3<- lm(polity_for1~death*agepre+polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m4<- lm(polity_for1~death*agepre+tenpre+polity_lag1,data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d1<-subset(d1,is.na(ccode)==FALSE &is.na(year)==FALSE)

d<-d1
m5<- lm(polity_for1~death*agepre+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE)
m6<- lm(polity_for1~death*agepre + polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))
a<-seq(quantile(d$agepre,0.1, na.rm=TRUE),quantile(d$agepre,0.9,na.rm=TRUE),1)
marg1_fe<-margins(m6,at=list(agepre=a),vcov=vc,variables="death")

d<-subset(nondem,is.na(polity_for1) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m7<- lm(polity_for1~death*agepre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m7,cluster=d$ccode,type="CR1")
se7<-sqrt(diag(vc))

d<-subset(d1,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m8<- lm(polity_for1~death*agepre+tenpre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m8,cluster=d$ccode,type="CR1")
se8<-sqrt(diag(vc))


# 5-Year DV (Appendix, Table 15)

nondem<-subset(data,polity_lag1<0)
nondem<-subset(nondem, is.na(death)==FALSE & is.na(agepre)==FALSE )

d5<-subset(nondem,death_for5==0 )
d5<-subset(d1,is.na(polity_for5)==FALSE)

d<-d5
m1<- lm(polity_for5~death*agepre,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d<-subset(d5,is.na(polity_lag1)==FALSE)
m2<- lm(polity_for5~death*agepre + polity_lag1,data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))
a<-seq(quantile(d$agepre,0.1, na.rm=TRUE),quantile(d$agepre,0.9,na.rm=TRUE),1)
marg5_nofe<-margins(m2,at=list(agepre=a),vcov=vc,variables="death")

d<-subset(nondem,is.na(polity_for5) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m3<- lm(polity_for5~death*agepre+polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))

d<-subset(d5,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m4<- lm(polity_for5~death*agepre+tenpre+polity_lag1,data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d5<-subset(d5,is.na(ccode)==FALSE &is.na(year)==FALSE)

d<-d5
m5<- lm(polity_for5~death*agepre+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))

d<-subset(d5,is.na(polity_lag1)==FALSE)
m6<- lm(polity_for5~death*agepre + polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))
a<-seq(quantile(d$agepre,0.1, na.rm=TRUE),quantile(d$agepre,0.9,na.rm=TRUE),1)
marg5_fe<-margins(m6,at=list(agepre=a),vcov=vc,variables="death")

d<-subset(nondem,is.na(polity_for5) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m7<- lm(polity_for5~death*agepre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m7,cluster=d$ccode,type="CR1")
se7<-sqrt(diag(vc))

d<-subset(d5,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m8<- lm(polity_for5~death*agepre+tenpre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m8,cluster=d$ccode,type="CR1")
se8<-sqrt(diag(vc))


# 10-Year DV (Appendix, Table 16)

nondem<-subset(data,polity_lag1<0)
nondem<-subset(nondem, is.na(death)==FALSE & is.na(agepre)==FALSE )

d10<-subset(nondem,death_for10==0 )
d10<-subset(d1,is.na(polity_for10)==FALSE)

d<-d10
m1<- lm(polity_for10~death*agepre,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d<-subset(d10,is.na(polity_lag1)==FALSE)
m2<- lm(polity_for10~death*agepre + polity_lag1,data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))
a<-seq(quantile(d$agepre,0.1, na.rm=TRUE),quantile(d$agepre,0.9,na.rm=TRUE),1)
marg10_nofe<-margins(m2,at=list(agepre=a),vcov=vc,variables="death")

d<-subset(nondem,is.na(polity_for10) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m3<- lm(polity_for10~death*agepre+polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))

d<-subset(d10,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m4<- lm(polity_for10~death*agepre+tenpre+polity_lag1,data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d10<-subset(d10,is.na(ccode)==FALSE &is.na(year)==FALSE)

d<-d10
m5<- lm(polity_for10~death*agepre+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se10<-sqrt(diag(vc))

d<-subset(d10,is.na(polity_lag1)==FALSE)
m6<- lm(polity_for10~death*agepre + polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))
a<-seq(quantile(d$agepre,0.1, na.rm=TRUE),quantile(d$agepre,0.9,na.rm=TRUE),1)
marg10_fe<-margins(m6,at=list(agepre=a),vcov=vc,variables="death")

d<-subset(nondem,is.na(polity_for10) ==FALSE & is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE)
m7<- lm(polity_for10~death*agepre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m7,cluster=d$ccode,type="CR1")
se7<-sqrt(diag(vc))

d<-subset(d10,is.na(polity_lag1)==FALSE & is.na(agepre)==FALSE &is.na(tenpre)==FALSE)
m8<- lm(polity_for10~death*agepre+tenpre+polity_lag1+factor(ccode)+factor(year),data=d)
vc<-vcovCR(m8,cluster=d$ccode,type="CR1")
se8<-sqrt(diag(vc))


### Plot for H1 : Figure 3 in main article 

confi<-function(x,length,lev){
  confi<-confint(subset(x, agepre==summary(x)$agepre[1]),level=lev)	
  for(i in c(summary(x)$agepre[1:length])){
    a<-i
    confi<-rbind(confi,confint(subset(x, agepre==a))	)
  }
  return(data.frame(confi))
}


pdf(file="figH1.pdf")
par(mfrow=c(3,2))
len<-nrow(summary(marg1_nofe))
confintty<-confi(marg1_nofe,len,0.9)
plot(x=summary(marg1_nofe)$agepre,y=as.numeric(summary(marg1_nofe)$AME),pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-2,1.5),cex=0.8, main="1-Year DV, No Fixed Effects")
segments(x0=summary(marg1_nofe)$agepre,
         y0=confintty$lower[1:len],
         x1=summary(marg1_nofe)$agepre,
         y1=confintty$upper[1:len])
abline(h=0)

len<-nrow(summary(marg1_fe))
confintty<-confi(marg1_fe,len,0.9)
plot(x=summary(marg1_fe)$agepre,y=as.numeric(summary(marg1_fe)$AME),pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-2,1.5),cex=0.8, main="1-Year DV, Fixed Effects")
segments(x0=summary(marg1_fe)$agepre,
         y0=confintty$lower[1:len],
         x1=summary(marg1_fe)$agepre,
         y1=confintty$upper[1:len])
abline(h=0)

len<-nrow(summary(marg5_nofe))
confintty<-confi(marg5_nofe,len,0.9)
plot(x=summary(marg5_nofe)$agepre,y=as.numeric(summary(marg5_nofe)$AME),pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-2,1.5),cex=0.8, main="5-Year DV, No Fixed Effects")
segments(x0=summary(marg5_nofe)$agepre,
         y0=confintty$lower[1:len],
         x1=summary(marg5_nofe)$agepre,
         y1=confintty$upper[1:len])
abline(h=0)

len<-nrow(summary(marg5_fe))
confintty<-confi(marg5_fe,len,0.9)
plot(x=summary(marg5_fe)$agepre,y=as.numeric(summary(marg5_fe)$AME),pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-2,1.5),cex=0.8, main="5-Year DV, Fixed Effects")
segments(x0=summary(marg5_fe)$agepre,
         y0=confintty$lower[1:len],
         x1=summary(marg5_fe)$agepre,
         y1=confintty$upper[1:len])
abline(h=0)

len<-nrow(summary(marg10_nofe))
confintty<-confi(marg10_nofe,len,0.9)
plot(x=summary(marg10_nofe)$agepre,y=as.numeric(summary(marg10_nofe)$AME),pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-2,1.5),cex=0.8, main="10-Year DV, No Fixed Effects")
segments(x0=summary(marg10_nofe)$agepre,
         y0=confintty$lower[1:len],
         x1=summary(marg10_nofe)$agepre,
         y1=confintty$upper[1:len])
abline(h=0)

len<-nrow(summary(marg10_fe))
confintty<-confi(marg10_fe,len,0.9)
plot(x=summary(marg10_fe)$agepre,y=as.numeric(summary(marg10_fe)$AME),pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-2,1.5),cex=0.8, main="10-Year DV, Fixed Effects")
segments(x0=summary(marg10_fe)$agepre,
         y0=confintty$lower[1:len],
         x1=summary(marg10_fe)$agepre,
         y1=confintty$upper[1:len])
abline(h=0)
dev.off()


######### Hypothesis 2a: Effect of leader death and age in personalist vs non-personalist regimes


# 1-Year DV (Appendix, Table 17)

## Basic Specification 

nondem<-subset(data,military_lag1==1 |monarchy_lag1==1|personal_lag1==1|party_lag1==1)
nondem<-subset(nondem, is.na(death)==FALSE & is.na(agepre)==FALSE & is.na(personal_lag1)==FALSE)

d1<-subset(nondem,death_for1==0 )
d1<-subset(d1,is.na(polity_for1)==FALSE & is.na(polity_lag1)==FALSE)
d1p<-subset(d1,personal_lag1==1)
d1n<-subset(d1,personal_lag1==0)

d<-d1p
m1<- lm(polity_for1~death*agepre + polity_lag1,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d1p<-subset(d1p,is.na(year)==FALSE)
d<-d1p
m2<- lm(polity_for1~death*agepre+ polity_lag1 +factor(year),data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))

d<-d1n
m3<- lm(polity_for1~death*agepre + polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))

d1n<-subset(d1n,is.na(year)==FALSE)
d<-d1n
m4<- lm(polity_for1~death*agepre+ polity_lag1+factor(year),data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-d1
m5<- lm(polity_for1~death*agepre*personal_lag1 + polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))
a<-seq(quantile(d$agepre,0.1, na.rm=TRUE),quantile(d$agepre,0.9,na.rm=TRUE),1)
marg1_nofe<-margins(m5,at=list(personal_lag1=0:1,agepre=a),vcov=vc,variables="death")


d1<-subset(d1,is.na(year)==FALSE)
d<-d1
m6<- lm(polity_for1~death*agepre*personal_lag1+ polity_lag1+factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))
marg1_fe<-margins(m6,at=list(personal_lag1=0:1,agepre=a),vcov=vc,variables="death")


## Tenure Specification (Appendix, Table 18)

nondem<-subset(data,military_lag1==1 |monarchy_lag1==1|personal_lag1==1|party_lag1==1)
nondem<-subset(nondem, is.na(death)==FALSE & is.na(agepre)==FALSE & is.na(personal_lag1)==FALSE)

d1<-subset(nondem,death_for1==0 )
d1<-subset(d1,is.na(polity_for1)==FALSE & is.na(polity_lag1)==FALSE & is.na(tenpre)==FALSE)
d1p<-subset(d1,personal_lag1==1)
d1n<-subset(d1,personal_lag1==0)

d<-d1p
m1<- lm(polity_for1~death*agepre + tenpre+ polity_lag1,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d1p<-subset(d1p,is.na(year)==FALSE)
d<-d1p
m2<- lm(polity_for1~death*agepre+tenpre+ polity_lag1+factor(year),data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))

d<-d1n
m3<- lm(polity_for1~death*agepre + tenpre+ polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))

d1n<-subset(d1n,is.na(year)==FALSE)
d<-d1n
m4<- lm(polity_for1~death*agepre+tenpre+polity_lag1+factor(year),data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-d1
m5<- lm(polity_for1~death*agepre*personal_lag1 +tenpre+ polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))

d1<-subset(d1,is.na(year)==FALSE)
d<-d1
m6<- lm(polity_for1~death*agepre*personal_lag1+tenpre+polity_lag1+ factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))

## Basic Specification, Full Sample (Appendix, Table 19)

nondem<-subset(data,military_lag1==1 |monarchy_lag1==1|personal_lag1==1|party_lag1==1)
nondem<-subset(nondem, is.na(death)==FALSE & is.na(agepre)==FALSE & is.na(personal_lag1)==FALSE)


d1<-subset(d1,is.na(polity_lag1)==FALSE)
d1p<-subset(d1,personal_lag1==1)
d1n<-subset(d1,personal_lag1==0)

d<-d1p
m1<- lm(polity_for1~death*agepre + polity_lag1,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d1p<-subset(d1p,is.na(year)==FALSE)
d<-d1p
m2<- lm(polity_for1~death*agepre+ polity_lag1 +factor(year),data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))

d<-d1n
m3<- lm(polity_for1~death*agepre + polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))

d1n<-subset(d1n,is.na(year)==FALSE)
d<-d1n
m4<- lm(polity_for1~death*agepre+ polity_lag1+factor(year),data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-d1
m5<- lm(polity_for1~death*agepre*personal_lag1 + polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))

d1<-subset(d1,is.na(year)==FALSE)
d<-d1
m6<- lm(polity_for1~death*agepre*personal_lag1+ polity_lag1+factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))


##### 5-Year DV

## Basic Specification (Appendix, Table 20)

nondem<-subset(data,military_lag1==1 |monarchy_lag1==1|personal_lag1==1|party_lag1==1)
nondem<-subset(nondem, is.na(death)==FALSE & is.na(agepre)==FALSE & is.na(personal_lag1)==FALSE)

d5<-subset(nondem,death_for5==0 )
d5<-subset(d5,is.na(polity_for5)==FALSE & is.na(polity_lag1)==FALSE)
d5p<-subset(d5,personal_lag1==1)
d5n<-subset(d5,personal_lag1==0)

d<-d5p
m1<- lm(polity_for5~death*agepre + polity_lag1,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d5p<-subset(d5p,is.na(year)==FALSE)
d<-d5p
m2<- lm(polity_for5~death*agepre+ polity_lag1 +factor(year),data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))

d<-d5n
m3<- lm(polity_for5~death*agepre + polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))

d5n<-subset(d5n,is.na(year)==FALSE)
d<-d5n
m4<- lm(polity_for5~death*agepre+ polity_lag1+factor(year),data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-d5
m5<- lm(polity_for5~death*agepre*personal_lag1 + polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))
marg5_nofe<-margins(m5,at=list(personal_lag1=0:1,agepre=a),vcov=vc,variables="death")


d5<-subset(d5,is.na(year)==FALSE)
d<-d5
m6<- lm(polity_for5~death*agepre*personal_lag1+ polity_lag1+factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))
marg5_fe<-margins(m6,at=list(personal_lag1=0:1,agepre=a),vcov=vc,variables="death")


## Tenure Specification (Appendix, Table 21)

nondem<-subset(data,military_lag1==1 |monarchy_lag1==1|personal_lag1==1|party_lag1==1)
nondem<-subset(nondem, is.na(death)==FALSE & is.na(agepre)==FALSE & is.na(personal_lag1)==FALSE)

d5<-subset(nondem,death_for5==0 )
d5<-subset(d5,is.na(polity_for5)==FALSE & is.na(polity_lag1)==FALSE & is.na(tenpre)==FALSE)
d5p<-subset(d5,personal_lag1==1)
d5n<-subset(d5,personal_lag1==0)

d<-d5p
m1<- lm(polity_for5~death*agepre + tenpre+ polity_lag1,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d5p<-subset(d5p,is.na(year)==FALSE)
d<-d5p
m2<- lm(polity_for5~death*agepre+tenpre+ polity_lag1+factor(year),data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))

d<-d5n
m3<- lm(polity_for5~death*agepre + tenpre+ polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))

d5n<-subset(d5n,is.na(year)==FALSE)
d<-d5n
m4<- lm(polity_for5~death*agepre+tenpre+polity_lag1+factor(year),data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-d5
m5<- lm(polity_for5~death*agepre*personal_lag1 +tenpre+ polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))

d5<-subset(d5,is.na(year)==FALSE)
d<-d5
m6<- lm(polity_for5~death*agepre*personal_lag1+tenpre+polity_lag1+ factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))


## Basic Specification, Full Sample  (Appendix, Table 22)

nondem<-subset(data,military_lag1==1 |monarchy_lag1==1|personal_lag1==1|party_lag1==1)
nondem<-subset(nondem, is.na(death)==FALSE & is.na(agepre)==FALSE & is.na(personal_lag1)==FALSE)


d5<-subset(d5,is.na(polity_lag1)==FALSE)
d5p<-subset(d5,personal_lag1==1)
d5n<-subset(d5,personal_lag1==0)

d<-d5p
m1<- lm(polity_for5~death*agepre + polity_lag1,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d5p<-subset(d5p,is.na(year)==FALSE)
d<-d5p
m2<- lm(polity_for5~death*agepre+ polity_lag1 +factor(year),data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))

d<-d5n
m3<- lm(polity_for5~death*agepre + polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))

d5n<-subset(d5n,is.na(year)==FALSE)
d<-d5n
m4<- lm(polity_for5~death*agepre+ polity_lag1+factor(year),data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-d5
m5<- lm(polity_for5~death*agepre*personal_lag1 + polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))

d5<-subset(d5,is.na(year)==FALSE)
d<-d5
m6<- lm(polity_for5~death*agepre*personal_lag1+ polity_lag1+factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))

### 10-Year DV

## Basic Specification (Appendix, Table 23)

nondem<-subset(data,military_lag1==1 |monarchy_lag1==1|personal_lag1==1|party_lag1==1)
nondem<-subset(nondem, is.na(death)==FALSE & is.na(agepre)==FALSE & is.na(personal_lag1)==FALSE)

d10<-subset(nondem,death_for10==0 )
d10<-subset(d10,is.na(polity_for10)==FALSE & is.na(polity_lag1)==FALSE)
d10p<-subset(d10,personal_lag1==1)
d10n<-subset(d10,personal_lag1==0)

d<-d10p
m1<- lm(polity_for10~death*agepre + polity_lag1,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d10p<-subset(d10p,is.na(year)==FALSE)
d<-d10p
m2<- lm(polity_for10~death*agepre+ polity_lag1 +factor(year),data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))

d<-d10n
m3<- lm(polity_for10~death*agepre + polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))

d10n<-subset(d10n,is.na(year)==FALSE)
d<-d10n
m4<- lm(polity_for10~death*agepre+ polity_lag1+factor(year),data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-d10
m5<- lm(polity_for10~death*agepre*personal_lag1 + polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))
marg10_nofe<-margins(m5,at=list(personal_lag1=0:1,agepre=a),vcov=vc,variables="death")


d10<-subset(d10,is.na(year)==FALSE)
d<-d10
m6<- lm(polity_for10~death*agepre*personal_lag1+ polity_lag1+factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))
marg10_fe<-margins(m6,at=list(personal_lag1=0:1,agepre=a),vcov=vc,variables="death")


## Tenure Specification (Appendix, Table 24)

nondem<-subset(data,military_lag1==1 |monarchy_lag1==1|personal_lag1==1|party_lag1==1)
nondem<-subset(nondem, is.na(death)==FALSE & is.na(agepre)==FALSE & is.na(personal_lag1)==FALSE)

d10<-subset(nondem,death_for10==0 )
d10<-subset(d10,is.na(polity_for10)==FALSE & is.na(polity_lag1)==FALSE & is.na(tenpre)==FALSE)
d10p<-subset(d10,personal_lag1==1)
d10n<-subset(d10,personal_lag1==0)

d<-d10p
m1<- lm(polity_for10~death*agepre + tenpre+ polity_lag1,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d10p<-subset(d10p,is.na(year)==FALSE)
d<-d10p
m2<- lm(polity_for10~death*agepre+tenpre+ polity_lag1+factor(year),data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))

d<-d10n
m3<- lm(polity_for10~death*agepre + tenpre+ polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))

d10n<-subset(d10n,is.na(year)==FALSE)
d<-d10n
m4<- lm(polity_for10~death*agepre+tenpre+polity_lag1+factor(year),data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-d10
m5<- lm(polity_for10~death*agepre*personal_lag1 +tenpre+ polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))

d10<-subset(d10,is.na(year)==FALSE)
d<-d10
m6<- lm(polity_for10~death*agepre*personal_lag1+tenpre+polity_lag1+ factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))

## Basic Specification, Full Sample (Appendix, Table 25)

nondem<-subset(data,military_lag1==1 |monarchy_lag1==1|personal_lag1==1|party_lag1==1)
nondem<-subset(nondem, is.na(death)==FALSE & is.na(agepre)==FALSE & is.na(personal_lag1)==FALSE)


d10<-subset(d10,is.na(polity_lag1)==FALSE)
d10p<-subset(d10,personal_lag1==1)
d10n<-subset(d10,personal_lag1==0)

d<-d10p
m1<- lm(polity_for10~death*agepre + polity_lag1,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

d10p<-subset(d10p,is.na(year)==FALSE)
d<-d10p
m2<- lm(polity_for10~death*agepre+ polity_lag1 +factor(year),data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))

d<-d10n
m3<- lm(polity_for10~death*agepre + polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))

d10n<-subset(d10n,is.na(year)==FALSE)
d<-d10n
m4<- lm(polity_for10~death*agepre+ polity_lag1+factor(year),data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

d<-d10
m5<- lm(polity_for10~death*agepre*personal_lag1 + polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))

d10<-subset(d10,is.na(year)==FALSE)
d<-d10
m6<- lm(polity_for10~death*agepre*personal_lag1+ polity_lag1+factor(year),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))


####FIGURE FOR H2a (Figure 4 in Main Article)

confi<-function(x,length,pers,lev){
	p<-pers
	confi<-confint(subset(x,personal_lag1==p & agepre==summary(x)$agepre[1]),level=lev)	
	for(i in c(summary(x)$agepre[1:length])){
		a<-i
		confi<-rbind(confi,confint(subset(x,personal_lag1==p & agepre==a))	)
	}
	return(data.frame(confi))
}


pdf(file="figH2a.pdf")
par(mfrow=c(3,2))

len<-nrow(summary(marg1_nofe))
l1<-len/2
l2<-l1+1
confintty0<-confi(marg1_nofe,l1,0,0.9)
confintty1<-confi(marg1_nofe,l1,1,0.9)
plot(x=summary(marg1_nofe)$agepre[1:l1],y=summary(marg1_nofe)$AME[1:l1],pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-10,10),cex=0.8, main="1-Year DV, No Fixed Effects")
segments(	x0=summary(marg1_nofe)$agepre[1:l1],
			y0=confintty0$lower[1:l1],
			x1=summary(marg1_nofe)$agepre[1:l1],
			y1=confintty0$upper[1:l1])
points(x=summary(marg1_nofe)$agepre[l2:len],y=summary(marg1_nofe)$AME[l2:len],pch=20,col="gray")
segments(	x0=summary(marg1_nofe)$agepre[1:l1],
			y0=confintty1$lower[1:l1],
			x1=summary(marg1_nofe)$agepre[1:l1],
			y1=confintty1$upper[1:l1],col="gray")
abline(h=0)
legend(41,9,c("Non-Personalist","Personalist"), pch=19,col=c("black","gray"),text.font=3,cex=0.8)



len<-nrow(summary(marg1_fe))
l1<-len/2
l2<-l1+1
confintty0<-confi(marg1_fe,l1,0,0.9)
confintty1<-confi(marg1_fe,l1,1,0.9)
plot(x=summary(marg1_fe)$agepre[1:l1],y=summary(marg1_fe)$AME[1:l1],pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-10,10),cex=0.8, main="1-Year DV, Fixed Effects")
segments(	x0=summary(marg1_nofe)$agepre[1:l1],
			y0=confintty0$lower[1:l1],
			x1=summary(marg1_nofe)$agepre[1:l1],
			y1=confintty0$upper[1:l1])
points(x=summary(marg1_fe)$agepre[l2:len],y=summary(marg1_fe)$AME[l2:len],pch=20,col="gray")
segments(	x0=summary(marg1_fe)$agepre[1:l1],
			y0=confintty1$lower[1:l1],
			x1=summary(marg1_fe)$agepre[1:l1],
			y1=confintty1$upper[1:l1],col="gray")
abline(h=0)
legend(41,9,c("Non-Personalist","Personalist"), pch=19,col=c("black","gray"),text.font=3,cex=0.8)


len<-nrow(summary(marg5_nofe))
l1<-len/2
l2<-l1+1
confintty0<-confi(marg5_nofe,l1,0,0.9)
confintty1<-confi(marg5_nofe,l1,1,0.9)
plot(x=summary(marg5_nofe)$agepre[1:l1],y=summary(marg5_nofe)$AME[1:l1],pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-10,10),cex=0.8, main="5-Year DV, No Fixed Effects")
segments(	x0=summary(marg5_nofe)$agepre[1:l1],
			y0=confintty0$lower[1:l1],
			x1=summary(marg5_nofe)$agepre[1:l1],
			y1=confintty0$upper[1:l1])
points(x=summary(marg5_nofe)$agepre[l2:len],y=summary(marg5_nofe)$AME[l2:len],pch=20,col="gray")
segments(	x0=summary(marg5_nofe)$agepre[1:l1],
			y0=confintty1$lower[1:l1],
			x1=summary(marg5_nofe)$agepre[1:l1],
			y1=confintty1$upper[1:l1],col="gray")
abline(h=0)
legend(41,9,c("Non-Personalist","Personalist"), pch=19,col=c("black","gray"),text.font=3,cex=0.8)


len<-nrow(summary(marg5_fe))
l1<-len/2
l2<-l1+1
confintty0<-confi(marg5_fe,l1,0,0.9)
confintty1<-confi(marg5_fe,l1,1,0.9)
plot(x=summary(marg5_fe)$agepre[1:l1],y=summary(marg5_fe)$AME[1:l1],pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-10,10),cex=0.8, main="5-Year DV, Fixed Effects")
segments(	x0=summary(marg5_fe)$agepre[1:l1],
			y0=confintty0$lower[1:l1],
			x1=summary(marg5_fe)$agepre[1:l1],
			y1=confintty0$upper[1:l1])
points(x=summary(marg5_fe)$agepre[l2:len],y=summary(marg5_fe)$AME[l2:len],pch=20,col="gray")
segments(	x0=summary(marg5_fe)$agepre[1:l1],
			y0=confintty1$lower[1:l1],
			x1=summary(marg5_fe)$agepre[1:l1],
			y1=confintty1$upper[1:l1],col="gray")
abline(h=0)
legend(41,9,c("Non-Personalist","Personalist"), pch=19,col=c("black","gray"),text.font=3,cex=0.8)



len<-nrow(summary(marg10_nofe))
l1<-len/2
l2<-l1+1
confintty0<-confi(marg10_nofe,l1,0,0.9)
confintty1<-confi(marg10_nofe,l1,1,0.9)
plot(x=summary(marg10_nofe)$agepre[1:l1],y=summary(marg10_nofe)$AME[1:l1],pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-10,10),cex=0.8, main="10-Year DV, No Fixed Effects")
segments(	x0=summary(marg10_nofe)$agepre[1:l1],
			y0=confintty0$lower[1:l1],
			x1=summary(marg10_nofe)$agepre[1:l1],
			y1=confintty0$upper[1:l1])
points(x=summary(marg10_nofe)$agepre[l2:len],y=summary(marg10_nofe)$AME[l2:len],pch=20,col="gray")
segments(	x0=summary(marg10_nofe)$agepre[1:l1],
			y0=confintty1$lower[1:l1],
			x1=summary(marg10_nofe)$agepre[1:l1],
			y1=confintty1$upper[1:l1],col="gray")
abline(h=0)
legend(41,9,c("Non-Personalist","Personalist"), pch=19,col=c("black","gray"),text.font=3,cex=0.8)


len<-nrow(summary(marg10_fe))
l1<-len/2
l2<-l1+1
confintty0<-confi(marg10_fe,l1,0,0.9)
confintty1<-confi(marg10_fe,l1,1,0.9)
plot(x=summary(marg10_fe)$agepre[1:l1],y=summary(marg10_fe)$AME[1:l1],pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-10,10),cex=0.8, main="10-Year DV, Fixed Effects")
segments(	x0=summary(marg10_fe)$agepre[1:l1],
			y0=confintty0$lower[1:l1],
			x1=summary(marg10_fe)$agepre[1:l1],
			y1=confintty0$upper[1:l1])
points(x=summary(marg10_fe)$agepre[l2:len],y=summary(marg10_fe)$AME[l2:len],pch=20,col="gray")
segments(	x0=summary(marg10_fe)$agepre[1:l1],
			y0=confintty1$lower[1:l1],
			x1=summary(marg10_fe)$agepre[1:l1],
			y1=confintty1$upper[1:l1],col="gray")
abline(h=0)
legend(41,9,c("Non-Personalist","Personalist"), pch=19,col=c("black","gray"),text.font=3,cex=0.8)

dev.off()

####  H2b Analyses: Effect of leader age on likelihood of removal


### Logit analyses (Appendix, Table 26)

nondem<-subset(data,military_lag1==1 |monarchy_lag1==1|personal_lag1==1|party_lag1==1)


per<-subset(nondem,personal_lag1==1)
mil<-subset(nondem,military_lag1==1)
mon<-subset(nondem,monarchy_lag1==1)
par<-subset(nondem,party_lag1==1)


d<-nondem
m1 <- glm(allcoup_bin~agepre+personal_lag1 + polity_lag1 + allcoup_count_lag10,data=nondem,family=binomial(link=logit))
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

m2 <- glm(allcoup_bin~agepre*personal_lag1 + polity_lag1 + allcoup_count_lag10,data=nondem,family=binomial(link=logit))
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))
plotty1<-m2

m3<-clogit(allcoup_bin~agepre*personal_lag1 + polity_lag1+allcoup_count_lag10+strata(year),data=nondem)
se3<-sqrt(diag(vcov(m3)))


m4 <- glm(otherirreg~agepre+personal_lag1  + polity_lag1+ allcoup_count_lag10,data=nondem,family=binomial(link=logit))
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))


m5 <- glm(otherirreg~agepre*personal_lag1  + polity_lag1+ allcoup_count_lag10 ,data=nondem,family=binomial(link=logit))
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))
plotty2<-m5

m6<-clogit(otherirreg~agepre*personal_lag1+ polity_lag1+ allcoup_count_lag10  +strata(year),data=nondem)
se6<-sqrt(diag(vcov(m6)))


##Logit analyses, tenure sepecification (Appendix, Table 27)

nondem<-subset(data,military_lag1==1 |monarchy_lag1==1|personal_lag1==1|party_lag1==1)


per<-subset(nondem,personal_lag1==1)
mil<-subset(nondem,military_lag1==1)
mon<-subset(nondem,monarchy_lag1==1)
par<-subset(nondem,party_lag1==1)


d<-nondem
m1 <- glm(allcoup_bin~agepre+personal_lag1 + tenpre+polity_lag1 + allcoup_count_lag10,data=nondem,family=binomial(link=logit))
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))

m2 <- glm(allcoup_bin~agepre*personal_lag1 + tenpre+polity_lag1 + allcoup_count_lag10,data=nondem,family=binomial(link=logit))
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))

m3<-clogit(allcoup_bin~agepre*personal_lag1 + tenpre+polity_lag1+allcoup_count_lag10+strata(year),data=nondem)
se3<-sqrt(diag(vcov(m3)))

m4 <- glm(otherirreg~agepre+personal_lag1  + tenpre+polity_lag1+allcoup_count_lag10,data=nondem,family=binomial(link=logit))
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))

m5 <- glm(otherirreg~agepre*personal_lag1  + tenpre+polity_lag1+allcoup_count_lag10 ,data=nondem,family=binomial(link=logit))
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))

m6<-clogit(otherirreg~agepre*personal_lag1+ tenpre+polity_lag1 +allcoup_count_lag10 +strata(year),data=nondem)
se6<-sqrt(diag(vcov(m3)))

## Figure for H2b (Figure 6 in Main Article)
pdf(file="figH2b.pdf")
par(mfrow=c(2,1))

p1<-mean(nondem$polity_lag1,na.rm=TRUE)
coup<-mean(nondem$allcoup_count_lag10,na.rm=TRUE)
nondemd<-subset(nondem,death==1)
age<-seq(quantile(nondemd$agepre,0.1, na.rm=TRUE),quantile(nondemd$agepre,0.9,na.rm=TRUE),1)

eff1<-predict(plotty1, newdata=data.frame(polity_lag1=p1,allcoup_count_lag10=coup,agepre=age,personal_lag1=1),type="link",se.fit=TRUE)
up <- eff1$fit + 2*eff1$se.fit
low <- eff1$fit - 2*eff1$se.fit
upper_p<- exp(up)/(1+exp(up))
lower_p<- exp(low)/(1+exp(low))
eff1_p<-predict(plotty1, newdata=data.frame(polity_lag1=p1,allcoup_count_lag10=coup,agepre=age,personal_lag1=1),type="response")

eff1<-predict(plotty1, newdata=data.frame(polity_lag1=p1,allcoup_count_lag10=coup,agepre=age,personal_lag1=0),type="link",se.fit=TRUE)
up <- eff1$fit + 2*eff1$se.fit
low <- eff1$fit - 2*eff1$se.fit
upper_n <- exp(up)/(1+exp(up))
lower_n<- exp(low)/(1+exp(low))
eff1_n<-predict(plotty1, newdata=data.frame(polity_lag1=p1,allcoup_count_lag10=coup,agepre=age,personal_lag1=0),type="response")


plot(age,eff1_p,ylim=c(0,0.4),type="l",ylab="Predicted Probability",xlab="Leader Age",main="DV: Attempted Coup")
cord.x<-c(age[1],age[1],age,age[length(age)],age[length(age)],rev(age))
cord.y<-c(upper_p[1],lower_p[1],lower_p,lower_p[length(age)],upper_p[length(age)],rev(upper_p))
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.6))
lines(age,eff1_n)
cord.x<-c(age[1],age[1],age,age[length(age)],age[length(age)],rev(age))
cord.y<-c(upper_n[1],lower_n[1],lower_n,lower_n[length(age)],upper_n[length(age)],rev(upper_n))
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))
legend(52,0.385,c("Personalist","Non-Personalist"), pch=19,col=c(rgb(0,0,0,alpha=0.6),rgb(0,0,0,alpha=0.3)),text.font=3,cex=0.8)


p1<-mean(nondem$polity_lag1,na.rm=TRUE)
coup<-mean(nondem$allcoup_count_lag10,na.rm=TRUE)
nondemd<-subset(nondem,death==1)
age<-seq(quantile(nondemd$agepre,0.1, na.rm=TRUE),quantile(nondemd$agepre,0.9,na.rm=TRUE),1)

eff1<-predict(plotty2, newdata=data.frame(polity_lag1=p1,allcoup_count_lag10=coup,agepre=age,personal_lag1=1),type="link",se.fit=TRUE)
up <- eff1$fit + 2*eff1$se.fit
low <- eff1$fit - 2*eff1$se.fit
upper_p<- exp(up)/(1+exp(up))
lower_p<- exp(low)/(1+exp(low))
eff1_p<-predict(plotty2, newdata=data.frame(polity_lag1=p1,allcoup_count_lag10=coup,agepre=age,personal_lag1=1),type="response")

eff1<-predict(plotty2, newdata=data.frame(polity_lag1=p1,allcoup_count_lag10=coup,agepre=age,personal_lag1=0),type="link",se.fit=TRUE)
up <- eff1$fit + 2*eff1$se.fit
low <- eff1$fit - 2*eff1$se.fit
upper_n <- exp(up)/(1+exp(up))
lower_n<- exp(low)/(1+exp(low))
eff1_n<-predict(plotty2, newdata=data.frame(polity_lag1=p1,allcoup_count_lag10=coup,agepre=age,personal_lag1=0),type="response")


plot(age,eff1_p,ylim=c(0,0.4),type="l",ylab="Predicted Probability",xlab="Leader Age",main="DV: Irregular Exit")
cord.x<-c(age[1],age[1],age,age[length(age)],age[length(age)],rev(age))
cord.y<-c(upper_p[1],lower_p[1],lower_p,lower_p[length(age)],upper_p[length(age)],rev(upper_p))
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.6))
lines(age,eff1_n)
cord.x<-c(age[1],age[1],age,age[length(age)],age[length(age)],rev(age))
cord.y<-c(upper_n[1],lower_n[1],lower_n,lower_n[length(age)],upper_n[length(age)],rev(upper_n))
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))
legend(52,0.385,c("Personalist","Non-Personalist"), pch=19,col=c(rgb(0,0,0,alpha=0.6),rgb(0,0,0,alpha=0.3)),text.font=3,cex=0.8)

dev.off()


#### H2b: Survival analysis 

datasurv<-read.csv("survdata_190123.csv")
datasurv<-datasurv[-c(1)]
datasurv$start<-datasurv$tenpre
datasurv$stop<-datasurv$start+1


## DV: Coup Attempts (Appendix, Table 28)
nondemsurv<-subset(datasurv,polity_lag1<0)

m1<-coxph(Surv(start,stop,allcoup_bin)~agepre+personal_lag1+polity_lag1,data=nondemsurv)
m2<-coxph(Surv(start,stop,allcoup_bin)~agepre*personal_lag1+polity_lag1,data=nondemsurv)
m3<-coxph(Surv(start,stop,allcoup_bin)~agepre*personal_lag1 +polity_lag1+ allcoup_count_lag10,data=nondemsurv)
m4<-coxph(Surv(start,stop,allcoup_bin)~agepre+personal_lag1+polity_lag1+factor(year),data=nondemsurv)
m5<-coxph(Surv(start,stop,allcoup_bin)~agepre*personal_lag1+polity_lag1+factor(year),data=nondemsurv)
m6<-coxph(Surv(start,stop,allcoup_bin)~agepre*personal_lag1 +polity_lag1+factor(year)+ allcoup_count_lag10,data=nondemsurv)

## DV: Coup Attempts (Appendix, Table 29)
m7<-coxph(Surv(start,stop,otherirreg)~agepre+personal_lag1+polity_lag1,data=nondemsurv)
m8<-coxph(Surv(start,stop,otherirreg)~agepre*personal_lag1+polity_lag1,data=nondemsurv)
m9<-coxph(Surv(start,stop,otherirreg)~agepre*personal_lag1+polity_lag1+allcoup_count_lag10,data=nondemsurv)
m10<-coxph(Surv(start,stop,otherirreg)~agepre+personal_lag1+polity_lag1+factor(year),data=nondemsurv)
m11<-coxph(Surv(start,stop,otherirreg)~agepre*personal_lag1+polity_lag1+factor(year),data=nondemsurv)
m12<-coxph(Surv(start,stop,otherirreg)~agepre*personal_lag1+polity_lag1+allcoup_count_lag10+factor(year),data=nondemsurv)



##H2b Survival Analysis Figure (Figure 7 in Appendix)

nondemsurvd<-subset(nondemsurv,death==1)	
al<-round(quantile(nondemsurvd$agepre,0.1, na.rm=TRUE))
ah<-round(quantile(nondemsurvd$agepre,0.9,na.rm=TRUE))

nondem<-subset(data,polity_lag1<0)	
nondemd<-subset(nondem,death==1)	
al<-round(quantile(nondemd$agepre,0.1, na.rm=TRUE))
ah<-round(quantile(nondemd$agepre,0.9,na.rm=TRUE))


pdf(file="figH2b_Surv.pdf")
par(mfrow=c(2,2))

xmax<-20

newdata_young<-data.frame(agepre=al, polity_lag1=mean(nondemsurv$polity_lag1,na.rm=TRUE),personal_lag1=0,allcoup_count_lag10=mean(nondemsurv$allcoup_count_lag10,na.rm=TRUE))
fit_young<-survfit(m3,newdata=newdata_young,conf.type="plain",conf.int=0.9)
newdata_old<-data.frame(agepre=ah,polity_lag1=mean(nondemsurv$polity_lag1,na.rm=TRUE), personal_lag1=0,allcoup_count_lag10=mean(nondemsurv$allcoup_count_lag10,na.rm=TRUE))
fit_old<-survfit(m3,newdata=newdata_old,conf.type="plain",conf.int=0.90)
plot(fit_young$time[1:xmax], fit_young$surv[1:xmax],lty=1,main = "DV: Coup Attempt",xlab="years in office",ylab="survival= no coup",type="l",ylim=c(0,1.2),xlim=c(fit_young$time[1],fit_young$time[xmax]),yaxt="n")
axis(2,at=c(0,0.2,0.4,0.6,0.8,1),labels=c("0","0.2","0.4","0.6","0.8","1"),las=2)
cord.x<-c(fit_young$time[1],fit_young$time[1],fit_young$time[1:xmax],fit_young$time[xmax],fit_young$time[xmax],rev(fit_young$time[1:xmax]))
cord.y<-c(fit_young$upper[1],fit_young$lower[1],fit_young$lower[1:xmax],fit_young$lower[xmax],fit_young$upper[xmax],rev(fit_young$upper[1:xmax]))
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.6))
lines(fit_old$time[1:xmax], fit_old$surv[1:xmax])
cord.x<-c(fit_old$time[1],fit_old$time[1],fit_old$time[1:xmax],fit_old$time[xmax],fit_old$time[xmax],rev(fit_old$time[1:xmax]))
cord.y<-c(fit_old$upper[1],fit_old$lower[1],fit_old$lower[1:xmax],fit_old$lower[xmax],fit_old$upper[xmax],rev(fit_old$upper[1:xmax]))
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))
abline(h=1.1)
legend(fit_old$time[xmax]-6,1.05,legend=c("Young","Old"),pch=19,col=c(rgb(0,0,0,alpha=0.6),rgb(0,0,0,alpha=0.3)))
text(3.5,1.17,c("Non-Personalist Regimes"),vfont=c("sans serif","bold italic"),pos=4)

newdata_young<-data.frame(agepre=al, polity_lag1=mean(nondemsurv$polity_lag1,na.rm=TRUE),personal_lag1=0,allcoup_count_lag10=mean(nondemsurv$allcoup_count_lag10,na.rm=TRUE))
fit_young<-survfit(m9,newdata=newdata_young,conf.type="plain",conf.int=0.9)
newdata_old<-data.frame(agepre=ah,polity_lag1=mean(nondemsurv$polity_lag1,na.rm=TRUE), personal_lag1=0,allcoup_count_lag10=mean(nondemsurv$allcoup_count_lag10,na.rm=TRUE))
fit_old<-survfit(m9,newdata=newdata_old,conf.type="plain",conf.int=0.90)
plot(fit_young$time[1:xmax], fit_young$surv[1:xmax],lty=1,main = "DV: Irregular Exit",xlab="years in office",ylab="survival= no irregular exit",type="l",ylim=c(0,1.2),xlim=c(fit_young$time[1],fit_young$time[xmax]),yaxt="n")
axis(2,at=c(0,0.2,0.4,0.6,0.8,1),labels=c("0","0.2","0.4","0.6","0.8","1"),las=2)
cord.x<-c(fit_young$time[1],fit_young$time[1],fit_young$time[1:xmax],fit_young$time[xmax],fit_young$time[xmax],rev(fit_young$time[1:xmax]))
cord.y<-c(fit_young$upper[1],fit_young$lower[1],fit_young$lower[1:xmax],fit_young$lower[xmax],fit_young$upper[xmax],rev(fit_young$upper[1:xmax]))
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.6))
lines(fit_old$time[1:xmax], fit_old$surv[1:xmax])
cord.x<-c(fit_old$time[1],fit_old$time[1],fit_old$time[1:xmax],fit_old$time[xmax],fit_old$time[xmax],rev(fit_old$time[1:xmax]))
cord.y<-c(fit_old$upper[1],fit_old$lower[1],fit_old$lower[1:xmax],fit_old$lower[xmax],fit_old$upper[xmax],rev(fit_old$upper[1:xmax]))
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))
abline(h=1.1)
legend(fit_old$time[xmax]-6,1.05,legend=c("Young","Old"),pch=19,col=c(rgb(0,0,0,alpha=0.6),rgb(0,0,0,alpha=0.3)))
text(3.5,1.17,c("Non-Personalist Regimes"),vfont=c("sans serif","bold italic"),pos=4)


newdata_young<-data.frame(agepre=al, polity_lag1=mean(nondemsurv$polity_lag1,na.rm=TRUE),personal_lag1=1,allcoup_count_lag10=mean(nondemsurv$allcoup_count_lag10,na.rm=TRUE))
fit_young<-survfit(m9,newdata=newdata_young,conf.type="plain",conf.int=0.9)
newdata_old<-data.frame(agepre=ah,polity_lag1=mean(nondemsurv$polity_lag1,na.rm=TRUE), personal_lag1=1,allcoup_count_lag10=mean(nondemsurv$allcoup_count_lag10,na.rm=TRUE))
fit_old<-survfit(m9,newdata=newdata_old,conf.type="plain",conf.int=0.90)
plot(fit_young$time[1:xmax], fit_young$surv[1:xmax],lty=1,main = "",xlab="years in office",ylab="survival= no coup",type="l",ylim=c(0,1.2),xlim=c(fit_young$time[1],fit_young$time[xmax]),yaxt="n")
axis(2,at=c(0,0.2,0.4,0.6,0.8,1),labels=c("0","0.2","0.4","0.6","0.8","1"),las=2)
cord.x<-c(fit_young$time[1],fit_young$time[1],fit_young$time[1:xmax],fit_young$time[xmax],fit_young$time[xmax],rev(fit_young$time[1:xmax]))
cord.y<-c(fit_young$upper[1],fit_young$lower[1],fit_young$lower[1:xmax],fit_young$lower[xmax],fit_young$upper[xmax],rev(fit_young$upper[1:xmax]))
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.6))
lines(fit_old$time[1:xmax], fit_old$surv[1:xmax])
cord.x<-c(fit_old$time[1],fit_old$time[1],fit_old$time[1:xmax],fit_old$time[xmax],fit_old$time[xmax],rev(fit_old$time[1:xmax]))
cord.y<-c(fit_old$upper[1],fit_old$lower[1],fit_old$lower[1:xmax],fit_old$lower[xmax],fit_old$upper[xmax],rev(fit_old$upper[1:xmax]))
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))
abline(h=1.1)
legend(fit_old$time[xmax]-6,1.05,legend=c("Young","Old"),pch=19,col=c(rgb(0,0,0,alpha=0.6),rgb(0,0,0,alpha=0.3)))
text(5,1.17,c("Personalist Regimes"),vfont=c("sans serif","bold italic"),pos=4)

newdata_young<-data.frame(agepre=al, polity_lag1=mean(nondemsurv$polity_lag1,na.rm=TRUE),personal_lag1=1,allcoup_count_lag10=mean(nondemsurv$allcoup_count_lag10,na.rm=TRUE))
fit_young<-survfit(m9,newdata=newdata_young,conf.type="plain",conf.int=0.9)
newdata_old<-data.frame(agepre=ah,polity_lag1=mean(nondemsurv$polity_lag1,na.rm=TRUE), personal_lag1=1,allcoup_count_lag10=mean(nondemsurv$allcoup_count_lag10,na.rm=TRUE))
fit_old<-survfit(m9,newdata=newdata_old,conf.type="plain",conf.int=0.90)
plot(fit_young$time[1:xmax], fit_young$surv[1:xmax],lty=1,main = "",xlab="years in office",ylab="survival= no irregular exit",type="l",ylim=c(0,1.2),xlim=c(fit_young$time[1],fit_young$time[xmax]),yaxt="n")
axis(2,at=c(0,0.2,0.4,0.6,0.8,1),labels=c("0","0.2","0.4","0.6","0.8","1"),las=2)
cord.x<-c(fit_young$time[1],fit_young$time[1],fit_young$time[1:xmax],fit_young$time[xmax],fit_young$time[xmax],rev(fit_young$time[1:xmax]))
cord.y<-c(fit_young$upper[1],fit_young$lower[1],fit_young$lower[1:xmax],fit_young$lower[xmax],fit_young$upper[xmax],rev(fit_young$upper[1:xmax]))
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.6))
lines(fit_old$time[1:xmax], fit_old$surv[1:xmax])
cord.x<-c(fit_old$time[1],fit_old$time[1],fit_old$time[1:xmax],fit_old$time[xmax],fit_old$time[xmax],rev(fit_old$time[1:xmax]))
cord.y<-c(fit_old$upper[1],fit_old$lower[1],fit_old$lower[1:xmax],fit_old$lower[xmax],fit_old$upper[xmax],rev(fit_old$upper[1:xmax]))
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))
abline(h=1.1)
legend(fit_old$time[xmax]-6,1.05,legend=c("Young","Old"),pch=19,col=c(rgb(0,0,0,alpha=0.6),rgb(0,0,0,alpha=0.3)))
text(1,1.17,c("Personalist Regimes"),vfont=c("sans serif","bold italic"),pos=4)

dev.off()


### Hypothesis 3: Effect of leader death and age, by level of economic development 


## Basic Specification (Appendix, Table 30)

nondem<-subset(data,polity_lag1<0)
nondem<-subset(nondem, is.na(death)==FALSE & is.na(agepre)==FALSE & is.na(gdppcln_lag1)==FALSE & is.na(polity_lag1)==FALSE)

gdp1<-quantile(nondem$gdppcln_lag1,0.1,na.rm=TRUE)
gdp2<-quantile(nondem$gdppcln_lag1,0.9,na.rm=TRUE)
a1<-quantile(nondem$agepre,0.1,na.rm=TRUE)
a2<-quantile(nondem$agepre,0.9,na.rm=TRUE)


d1<-subset(nondem,death_for1==0 & is.na(polity_for1)==FALSE )

d<-d1
m1<- lm(polity_for1~death*agepre*gdppcln_lag1 + polity_lag1,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))
a<-seq(quantile(d$agepre,0.1, na.rm=TRUE),quantile(d$agepre,0.9,na.rm=TRUE),1)
marg1_nofe<-margins(m1,at=list(gdppcln_lag1=c(gdp1,gdp2),agepre=a),vcov=vc,variables="death")

d1<-subset(d1,is.na(year)==FALSE & is.na(ccode)==FALSE)
d<-d1
m2<- lm(polity_for1~death*agepre*gdppcln_lag1+ polity_lag1+factor(year)+factor(ccode),data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))
marg1_fe<-margins(m2,at=list(gdppcln_lag1=c(gdp1,gdp2),agepre=a),vcov=vc,variables="death")


d5<-subset(nondem,death_for5==0 & is.na(polity_for5)==FALSE )

d<-d5
m3<- lm(polity_for5~death*agepre*gdppcln_lag1 + polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))
a<-seq(quantile(d$agepre,0.1, na.rm=TRUE),quantile(d$agepre,0.9,na.rm=TRUE),1)
marg5_nofe<-margins(m3,at=list(gdppcln_lag1=c(gdp1,gdp2),agepre=a),vcov=vc,variables="death")


d5<-subset(d5,is.na(year)==FALSE & is.na(ccode)==FALSE)
d<-d5
m4<- lm(polity_for5~death*agepre*gdppcln_lag1+ polity_lag1+factor(year)+factor(ccode),data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))
marg5_fe<-margins(m4,at=list(gdppcln_lag1=c(gdp1,gdp2),agepre=a),vcov=vc,variables="death")


d10<-subset(nondem,death_for10==0 & is.na(polity_for10)==FALSE )

d<-d10
m5<- lm(polity_for10~death*agepre*gdppcln_lag1 + polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))
a<-seq(quantile(d$agepre,0.1, na.rm=TRUE),quantile(d$agepre,0.9,na.rm=TRUE),1)
marg10_nofe<-margins(m5,at=list(gdppcln_lag1=c(gdp1,gdp2),agepre=a),vcov=vc,variables="death")


d10<-subset(d10,is.na(year)==FALSE & is.na(ccode)==FALSE)
d<-d10
m6<- lm(polity_for10~death*agepre*gdppcln_lag1+ polity_lag1+factor(year)+factor(ccode),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))
marg10_fe<-margins(m6,at=list(gdppcln_lag1=c(gdp1,gdp2),agepre=a),vcov=vc,variables="death")


## Tenure Specification  (Appendix, Table 31)
 

nondem<-subset(data,polity_lag1<0)
nondem<-subset(nondem, is.na(death)==FALSE & is.na(agepre)==FALSE & is.na(gdppcln_lag1)==FALSE & is.na(polity_lag1)==FALSE & is.na(tenpre)==FALSE)

gdp1<-quantile(nondem$gdppcln_lag1,0.1,na.rm=TRUE)
gdp2<-quantile(nondem$gdppcln_lag1,0.9,na.rm=TRUE)

d1<-subset(nondem,death_for1==0 & is.na(polity_for1)==FALSE )

d<-d1
m1<- lm(polity_for1~death*agepre*gdppcln_lag1 + tenpre+polity_lag1,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))


d1<-subset(d1,is.na(year)==FALSE & is.na(ccode)==FALSE)
d<-d1
m2<- lm(polity_for1~death*agepre*gdppcln_lag1+  tenpre+polity_lag1+factor(year)+factor(ccode),data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))


d5<-subset(nondem,death_for5==0 & is.na(polity_for5)==FALSE )

d<-d5
m3<- lm(polity_for5~death*agepre*gdppcln_lag1 +  tenpre+polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))

d5<-subset(d5,is.na(year)==FALSE & is.na(ccode)==FALSE)
d<-d5
m4<- lm(polity_for5~death*agepre*gdppcln_lag1+  tenpre+polity_lag1+factor(year)+factor(ccode),data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))


d10<-subset(nondem,death_for10==0 & is.na(polity_for10)==FALSE )

d<-d10
m5<- lm(polity_for10~death*agepre*gdppcln_lag1 +  tenpre+polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))


d10<-subset(d10,is.na(year)==FALSE & is.na(ccode)==FALSE)
d<-d10
m6<- lm(polity_for10~death*agepre*gdppcln_lag1+  tenpre+polity_lag1+factor(year)+factor(ccode),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))


## Full Sample (Appendix, Table 32)

nondem<-subset(data,polity_lag1<0)
nondem<-subset(nondem, is.na(death)==FALSE & is.na(agepre)==FALSE & is.na(gdppcln_lag1)==FALSE & is.na(polity_lag1)==FALSE)

gdp1<-quantile(nondem$gdppcln_lag1,0.1,na.rm=TRUE)
gdp2<-quantile(nondem$gdppcln_lag1,0.9,na.rm=TRUE)

d1<-subset(nondem,is.na(polity_for1)==FALSE )

d<-d1
m1<- lm(polity_for1~death*agepre*gdppcln_lag1 + polity_lag1,data=d)
vc<-vcovCR(m1,cluster=d$ccode,type="CR1")
se1<-sqrt(diag(vc))


d1<-subset(d1,is.na(year)==FALSE & is.na(ccode)==FALSE)
d<-d1
m2<- lm(polity_for1~death*agepre*gdppcln_lag1+ polity_lag1+factor(year)+factor(ccode),data=d)
vc<-vcovCR(m2,cluster=d$ccode,type="CR1")
se2<-sqrt(diag(vc))


d5<-subset(nondem, is.na(polity_for5)==FALSE )

d<-d5
m3<- lm(polity_for5~death*agepre*gdppcln_lag1 + polity_lag1,data=d)
vc<-vcovCR(m3,cluster=d$ccode,type="CR1")
se3<-sqrt(diag(vc))


d5<-subset(d5,is.na(year)==FALSE & is.na(ccode)==FALSE)
d<-d5
m4<- lm(polity_for5~death*agepre*gdppcln_lag1+ polity_lag1+factor(year)+factor(ccode),data=d)
vc<-vcovCR(m4,cluster=d$ccode,type="CR1")
se4<-sqrt(diag(vc))


d10<-subset(nondem,is.na(polity_for10)==FALSE )

d<-d10
m5<- lm(polity_for10~death*agepre*gdppcln_lag1 + polity_lag1,data=d)
vc<-vcovCR(m5,cluster=d$ccode,type="CR1")
se5<-sqrt(diag(vc))


d10<-subset(d10,is.na(year)==FALSE & is.na(ccode)==FALSE)
d<-d10
m6<- lm(polity_for10~death*agepre*gdppcln_lag1+ polity_lag1+factor(year)+factor(ccode),data=d)
vc<-vcovCR(m6,cluster=d$ccode,type="CR1")
se6<-sqrt(diag(vc))



## Plot for H3 (Figure 8 in Main Article)


confi<-function(x,length,gdp,lev){
	g<-gdp
	confi<-confint(subset(x,gdppcln_lag1==g & agepre==summary(x)$agepre[1]),level=lev)	
	for(i in c(summary(x)$agepre[1:length])){
		a<-i
		confi<-rbind(confi,confint(subset(x,gdppcln_lag1==g & agepre==a))	)
	}
	return(data.frame(confi))
}


pdf(file="figH3.pdf")
par(mfrow=c(3,2))

len<-nrow(summary(marg1_nofe))
l1<-len/2
l2<-l1+1
confintty0<-confi(marg1_nofe,l1,gdp1,0.9)
confintty1<-confi(marg1_nofe,l1,gdp2,0.9)
plot(x=summary(marg1_nofe)$agepre[1:l1],y=summary(marg1_nofe)$AME[1:l1],pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-10,10),cex=0.8, main="1-Year DV, No Fixed Effects")
segments(	x0=summary(marg1_nofe)$agepre[1:l1],
			y0=confintty0$lower[1:l1],
			x1=summary(marg1_nofe)$agepre[1:l1],
			y1=confintty0$upper[1:l1])
points(x=summary(marg1_nofe)$agepre[l2:len],y=summary(marg1_nofe)$AME[l2:len],pch=20,col="gray")
segments(	x0=summary(marg1_nofe)$agepre[1:l1],
			y0=confintty1$lower[1:l1],
			x1=summary(marg1_nofe)$agepre[1:l1],
			y1=confintty1$upper[1:l1],col="gray")
abline(h=0)
legend(60,9,c("Low-Development","High-Development"), pch=19,col=c("black","gray"),text.font=3,cex=0.8)

len<-nrow(summary(marg1_fe))
l1<-len/2
l2<-l1+1
confintty0<-confi(marg1_fe,l1,gdp1,0.9)
confintty1<-confi(marg1_fe,l1,gdp2,0.9)
plot(x=summary(marg1_fe)$agepre[1:l1],y=summary(marg1_fe)$AME[1:l1],pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-10,10),cex=0.8, main="1-Year DV, Fixed Effects")
segments(	x0=summary(marg1_fe)$agepre[1:l1],
			y0=confintty0$lower[1:l1],
			x1=summary(marg1_fe)$agepre[1:l1],
			y1=confintty0$upper[1:l1])
points(x=summary(marg1_fe)$agepre[l2:len],y=summary(marg1_fe)$AME[l2:len],pch=20,col="gray")
segments(	x0=summary(marg1_fe)$agepre[1:l1],
			y0=confintty1$lower[1:l1],
			x1=summary(marg1_fe)$agepre[1:l1],
			y1=confintty1$upper[1:l1],col="gray")
abline(h=0)
legend(60,9,c("Low-Development","High-Development"), pch=19,col=c("black","gray"),text.font=3,cex=0.8)

len<-nrow(summary(marg5_nofe))
l1<-len/2
l2<-l1+1
confintty0<-confi(marg5_nofe,l1,gdp1,0.9)
confintty1<-confi(marg5_nofe,l1,gdp2,0.9)
plot(x=summary(marg5_nofe)$agepre[1:l1],y=summary(marg5_nofe)$AME[1:l1],pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-10,10),cex=0.8, main="5-Year DV, No Fixed Effects")
segments(	x0=summary(marg5_nofe)$agepre[1:l1],
			y0=confintty0$lower[1:l1],
			x1=summary(marg5_nofe)$agepre[1:l1],
			y1=confintty0$upper[1:l1])
points(x=summary(marg5_nofe)$agepre[l2:len],y=summary(marg5_nofe)$AME[l2:len],pch=20,col="gray")
segments(	x0=summary(marg5_nofe)$agepre[1:l1],
			y0=confintty1$lower[1:l1],
			x1=summary(marg5_nofe)$agepre[1:l1],
			y1=confintty1$upper[1:l1],col="gray")
abline(h=0)
legend(59.4,9,c("Low-Development","High-Development"), pch=19,col=c("black","gray"),text.font=3,cex=0.8)

len<-nrow(summary(marg5_fe))
l1<-len/2
l2<-l1+1
confintty0<-confi(marg5_fe,l1,gdp1,0.9)
confintty1<-confi(marg5_fe,l1,gdp2,0.9)
plot(x=summary(marg5_fe)$agepre[1:l1],y=summary(marg5_fe)$AME[1:l1],pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-10,10),cex=0.8, main="5-Year DV, Fixed Effects")
segments(	x0=summary(marg5_fe)$agepre[1:l1],
			y0=confintty0$lower[1:l1],
			x1=summary(marg5_fe)$agepre[1:l1],
			y1=confintty0$upper[1:l1])
points(x=summary(marg5_fe)$agepre[l2:len],y=summary(marg5_fe)$AME[l2:len],pch=20,col="gray")
segments(	x0=summary(marg5_fe)$agepre[1:l1],
			y0=confintty1$lower[1:l1],
			x1=summary(marg5_fe)$agepre[1:l1],
			y1=confintty1$upper[1:l1],col="gray")
abline(h=0)
legend(59.4,9,c("Low-Development","High-Development"), pch=19,col=c("black","gray"),text.font=3,cex=0.8)

len<-nrow(summary(marg10_nofe))
l1<-len/2
l2<-l1+1
confintty0<-confi(marg10_nofe,l1,gdp1,0.9)
confintty1<-confi(marg10_nofe,l1,gdp2,0.9)
plot(x=summary(marg10_nofe)$agepre[1:l1],y=summary(marg10_nofe)$AME[1:l1],pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-10,10),cex=0.8, main="10-Year DV, No Fixed Effects")
segments(	x0=summary(marg10_nofe)$agepre[1:l1],
			y0=confintty0$lower[1:l1],
			x1=summary(marg10_nofe)$agepre[1:l1],
			y1=confintty0$upper[1:l1])
points(x=summary(marg10_nofe)$agepre[l2:len],y=summary(marg10_nofe)$AME[l2:len],pch=20,col="gray")
segments(	x0=summary(marg10_nofe)$agepre[1:l1],
			y0=confintty1$lower[1:l1],
			x1=summary(marg10_nofe)$agepre[1:l1],
			y1=confintty1$upper[1:l1],col="gray")
abline(h=0)
legend(58.4,9,c("Low-Development","High-Development"), pch=19,col=c("black","gray"),text.font=3,cex=0.8)

len<-nrow(summary(marg10_fe))
l1<-len/2
l2<-l1+1
confintty0<-confi(marg10_fe,l1,gdp1,0.9)
confintty1<-confi(marg10_fe,l1,gdp2,0.9)
plot(x=summary(marg10_fe)$agepre[1:l1],y=summary(marg10_fe)$AME[1:l1],pch=20,xlab="Age",ylab="Average Marginal Effect of Death",ylim=c(-10,10),cex=0.8, main="10-Year DV, Fixed Effects")
segments(	x0=summary(marg10_fe)$agepre[1:l1],
			y0=confintty0$lower[1:l1],
			x1=summary(marg10_fe)$agepre[1:l1],
			y1=confintty0$upper[1:l1])
points(x=summary(marg10_fe)$agepre[l2:len],y=summary(marg10_fe)$AME[l2:len],pch=20,col="gray")
segments(	x0=summary(marg10_fe)$agepre[1:l1],
			y0=confintty1$lower[1:l1],
			x1=summary(marg10_fe)$agepre[1:l1],
			y1=confintty1$upper[1:l1],col="gray")
abline(h=0)
legend(58.4,9,c("Low-Development","High-Development"), pch=19,col=c("black","gray"),text.font=3,cex=0.8)

dev.off()





####Descriptive Plots

nondem<-subset(data,democracy_lag1==0)

p1<-subset(nondem,death_for1==0)
p5<-subset(nondem,death_for5==0)
p10<-subset(nondem,death_for10==0)

p1d<-subset(p1,death==1)
p1nd<-subset(p1,death==0)



per<-subset(p1d,personal_lag1==1)
per<-per[order(per$agepre),]

q<-quantile(p1$gdppcln_lag1,.75,na.rm=TRUE) 
high<-subset(p1d,gdppcln_lag1>q)
high<-high[order(high$agepre),]

q<-quantile(p1$gdppcln_lag1,.25,na.rm=TRUE) 
low<-subset(p1d,gdppcln_lag1<q)
low<-low[order(low$agepre),]

c<-per$ccode
y<-per$year
n<-per$ccname
a<-per$agepre

## Personalist regimes (Appendix, Figure 1)

pdf(file="Personalist.pdf",8,8)
par(mfrow=c(3,3))
subby<-subset(data,ccode==c[1])
subby<-subset(subby,year>= y[1] - 10)
subby<-subset(subby,year<= y[1] + 10)
n[1]
a[1]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[1]-10,y[1]+10),type="l",ylab="Polity",xlab="Year",main="Mauritania: Age = 43")
abline(v=y[1],lty=1,col="red")	
dem<-subset(subby,democracy==1)

subby<-subset(data,ccode==c[2])
subby<-subset(subby,year>= y[2] - 10)
subby<-subset(subby,year<= y[2] + 10)
n[2]
a[2]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[2]-10,y[2]+10),type="l",ylab="Polity",xlab="Year",main="Iraq: Age = 45")
abline(v=y[2],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[3])
subby<-subset(subby,year>= y[3] - 10)
subby<-subset(subby,year<= y[3] + 10)
n[3]
a[3]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[3]-10,y[3]+10),type="l",ylab="Polity",xlab="Year",main="Bolivia: Age = 50")
abline(v=y[3],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
abline(v=dem$year[1],lty=1,col="black")

subby<-subset(data,ccode==c[4])
subby<-subset(subby,year>= y[4] - 10)
subby<-subset(subby,year<= y[4] + 10)
n[4]
a[4]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[4]-10,y[4]+10),type="l",ylab="Polity",xlab="Year",main="Nicaragua: Age = 57")
abline(v=y[4],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[5])
subby<-subset(subby,year>= y[5] - 10)
subby<-subset(subby,year<= y[5] + 10)
n[5]
a[5]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[5]-10,y[5]+10),type="l",ylab="Polity",xlab="Year",main="Haiti: Age = 64")
abline(v=y[5],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[6])
subby<-subset(subby,year>= y[6] - 10)
subby<-subset(subby,year<= y[6] + 10)
n[6]
a[6]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[6]-10,y[6]+10),type="l",ylab="Polity",xlab="Year",main="Portugal: Age = 79")
abline(v=y[6],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
cord.x<-c(dem$year[1],dem$year[1],dem$year[3],dem$year[3])
cord.y<-c(-11,11,11,-11)
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))


subby<-subset(data,ccode==c[7])
subby<-subset(subby,year>= y[7] - 10)
subby<-subset(subby,year<= y[7] + 10)
n[7]
a[7]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[7]-10,y[7]+10),type="l",ylab="Polity",xlab="Year",main="Spain: Age=83")
abline(v=y[7],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
cord.x<-c(dem$year[1],dem$year[1],dem$year[9],dem$year[9])
cord.y<-c(-11,11,11,-11)
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))

dev.off()

nondem<-subset(data,democracy_lag1==0)

p1<-subset(nondem,death_for1==0)
p5<-subset(nondem,death_for5==0)
p10<-subset(nondem,death_for10==0)

p1d<-subset(p1,death==1)
p1nd<-subset(p1,death==0)

mon<-subset(p1d,monarchy_lag1==1)
mon<-mon[order(mon$agepre),]
mil<-subset(p1d,military_lag1==1)
mil<-mil[order(mil$agepre),]
par<-subset(p1d,party_lag1==1)
par<-par[order(par$agepre),]


### Monarchy (Appendix, Figure 2)

c<-mon$ccode
y<-mon$year
n<-mon$ccname
a<-mon$agepre

pdf(file="Monarchy.pdf",6,6)
par(mfrow=c(4,3))

subby<-subset(data,ccode==c[2])
subby<-subset(subby,year>= y[2] - 10)
subby<-subset(subby,year<= y[2] + 10)
n[2]
a[2]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[2]-10,y[2]+10),type="l",ylab="Polity",xlab="Year",main="Morocco: Age = 52")
abline(v=y[2],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[3])
subby<-subset(subby,year>= y[3] - 10)
subby<-subset(subby,year<= y[3] + 10)
n[3]
a[3]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[3]-10,y[3]+10),type="l",ylab="Polity",xlab="Year",main="Nepal: Age = 52")
abline(v=y[3],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[4])
subby<-subset(subby,year>= y[4] - 10)
subby<-subset(subby,year<= y[4] + 10)
n[4]
a[4]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[4]-10,y[4]+10),type="l",ylab="Polity",xlab="Year",main="Jordan: Age = 64")
abline(v=y[4],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[5])
subby<-subset(subby,year>= y[5] - 10)
subby<-subset(subby,year<= y[5] + 10)
n[5]
a[5]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[5]-10,y[5]+10),type="l",ylab="Polity",xlab="Year",main="Kuwait: Age = 64")
abline(v=y[5],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[1])
subby<-subset(subby,year>= y[1] - 10)
subby<-subset(subby,year<= y[1] + 10)
n[1]
a[1]
n[6]
a[6]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[6]-10,y[1]+10),type="l",ylab="Polity",xlab="Year",main="Nepal: Age = 66, 49")
abline(v=y[1],lty=1,col="red")	
abline(v=y[6],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
subby<-subset(data,ccode==c[6])
subby<-subset(subby,year>= y[6] - 10)
subby<-subset(subby,year<= y[6] + 10)
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[7])
subby<-subset(subby,year>= y[7] - 10)
subby<-subset(subby,year<= y[7] + 10)
n[7]
a[7]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[7]-10,y[7]+10),type="l",ylab="Polity",xlab="Year",main="Saudi Arabia: Age = 69")
abline(v=y[7],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[8])
subby<-subset(subby,year>= y[8] - 10)
subby<-subset(subby,year<= y[8] + 10)
n[8]
a[8]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[8]-10,y[8]+10),type="l",ylab="Polity",xlab="Year",main="Morocco: Age = 70")
abline(v=y[8],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[9])
subby<-subset(subby,year>= y[9] - 10)
subby<-subset(subby,year<= y[9] + 10)
n[9]
a[9]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[9]-10,y[9]+10),type="l",ylab="Polity",xlab="Year",main="Kuwait: Age = 70")
abline(v=y[9],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)


subby<-subset(data,ccode==c[10])
subby<-subset(subby,year>= y[10] - 10)
subby<-subset(subby,year<= y[10] + 10)
n[10]
a[10]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[10]-10,y[10]+10),type="l",ylab="Polity",xlab="Year",main="Saudi Arabia: Age = 71")
abline(v=y[10],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[11])
subby<-subset(subby,year>= y[11] - 10)
subby<-subset(subby,year<= y[11] + 10)
n[11]
a[11]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[11]-10,y[11]+10),type="l",ylab="Polity",xlab="Year",main="Yemen: Age = 71")
abline(v=y[11],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[12])
subby<-subset(subby,year>= y[12] - 10)
subby<-subset(subby,year<= y[12] + 10)
n[12]
a[12]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[12]-10,y[12]+10),type="l",ylab="Polity",xlab="Year",main="Swaziland: Age = 83")
abline(v=y[12],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

dev.off()


## Military Regimes (Appendix, Figure 3)


c<-mil$ccode
y<-mil$year
n<-mil$ccname
a<-mil$agepre

pdf(file="Military.pdf",5,5)
par(mfrow=c(2,2))

subby<-subset(data,ccode==c[1])
subby<-subset(subby,year>= y[1] - 10)
subby<-subset(subby,year<= y[1] + 10)
n[1]
a[1]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[1]-10,y[1]+10),type="l",ylab="Polity",xlab="Year",main="Panama: Age = 52")
abline(v=y[1],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
cord.x<-c(dem$year[1],dem$year[1],dem$year[3],dem$year[3])
cord.y<-c(-11,11,11,-11)
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))

subby<-subset(data,ccode==c[2])
subby<-subset(subby,year>= y[2] - 10)
subby<-subset(subby,year<= y[2] + 10)
n[2]
a[2]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[2]-10,y[2]+10),type="l",ylab="Polity",xlab="Year",main="Nigeria: Age = 55")
abline(v=y[2],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
cord.x<-c(dem$year[1],dem$year[1],dem$year[10],dem$year[10])
cord.y<-c(-11,11,11,-11)
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))

subby<-subset(data,ccode==c[3])
subby<-subset(subby,year>= y[3] - 10)
subby<-subset(subby,year<= y[3] + 10)
n[3]
a[3]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[3]-10,y[3]+10),type="l",ylab="Polity",xlab="Year",main="Thailand: Age = 55")
abline(v=y[3],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[4])
subby<-subset(subby,year>= y[4] - 10)
subby<-subset(subby,year<= y[4] + 10)
n[4]
a[4]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[4]-10,y[4]+10),type="l",ylab="Polity",xlab="Year",main="Niger: Age = 56")
abline(v=y[4],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
cord.x<-c(dem$year[1],dem$year[1],dem$year[3],dem$year[3])
cord.y<-c(-11,11,11,-11)
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))

dev.off()


## Party Regimes (Appendix, Figure 4)


c<-par$ccode
y<-par$year
n<-par$ccname
a<-par$agepre

pdf(file="Party.pdf",12,15)
par(mfrow=c(8,4))

subby<-subset(data,ccode==c[1])
subby<-subset(subby,year>= y[1] - 10)
subby<-subset(subby,year<= y[1] + 10)
n[1]
a[1]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[1]-10,y[1]+10),type="l",ylab="Polity",xlab="Year",main="Algeria: Age = 51")
abline(v=y[1],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[2])
subby<-subset(subby,year>= y[2] - 10)
subby<-subset(subby,year<= y[2] + 10)
n[2]
a[2]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[2]-10,y[2]+10),type="l",ylab="Polity",xlab="Year",main="Egypt: Age = 52")
abline(v=y[2],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[3])
subby<-subset(subby,year>= y[3] - 10)
subby<-subset(subby,year<= y[3] + 10)
n[3]
a[3]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[3]-10,y[3]+10),type="l",ylab="Polity",xlab="Year",main="Mozambique: Age = 53")
abline(v=y[3],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[4])
subby<-subset(subby,year>= y[4] - 10)
subby<-subset(subby,year<= y[4] + 10)
n[4]
a[4]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[4]-10,y[4]+10),type="l",ylab="Polity",xlab="Year",main="Malaysia: Age = 54")
abline(v=y[4],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[5])
subby<-subset(subby,year>= y[5] - 10)
subby<-subset(subby,year<= y[16] + 10)
n[5]
a[5]
n[16]
a[16]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[5]-10,y[16]+10),type="l",ylab="Polity",xlab="Year",main="Czechoslovakia: Age = 57, 73")
abline(v=y[5],lty=1,col="red")	
abline(v=y[16],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
subby<-subset(data,ccode==c[16])
subby<-subset(subby,year>= y[16] - 10)
subby<-subset(subby,year<= y[16] + 10)
dem<-subset(subby,democracy==1)
nrow(dem)


subby<-subset(data,ccode==c[6])
subby<-subset(subby,year>= y[6] - 10)
subby<-subset(subby,year<= y[6] + 10)
n[6]
a[6]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[6]-10,y[6]+10),type="l",ylab="Polity",xlab="Year",main="Angola: Age = 57")
abline(v=y[6],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[7])
subby<-subset(subby,year>= y[7] - 10)
subby<-subset(subby,year<= y[7] + 10)
n[7]
a[7]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[7]-10,y[7]+10),type="l",ylab="Polity",xlab="Year",main="Mongolia: Age = 57")
abline(v=y[7],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[8])
subby<-subset(subby,year>= y[8] - 10)
subby<-subset(subby,year<= y[8] + 10)
n[8]
a[8]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[8]-10,y[8]+10),type="l",ylab="Polity",xlab="Year",main="Botswana: Age = 59")
abline(v=y[8],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[9])
subby<-subset(subby,year>= y[9] - 10)
subby<-subset(subby,year<= y[9] + 10)
n[9]
a[9]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[9]-10,y[9]+10),type="l",ylab="Polity",xlab="Year",main="Guinea: Age = 62")
abline(v=y[9],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[10])
subby<-subset(subby,year>= y[10] - 10)
subby<-subset(subby,year<= y[10] + 10)
n[10]
a[10]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[10]-10,y[10]+10),type="l",ylab="Polity",xlab="Year",main="Poland: Age = 64")
abline(v=y[10],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[11])
subby<-subset(subby,year>= y[11] - 10)
subby<-subset(subby,year<= y[11] + 10)
n[11]
a[11]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[11]-10,y[11]+10),type="l",ylab="Polity",xlab="Year",main="Romania: Age = 64")
abline(v=y[11],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[12])
subby<-subset(subby,year>= y[12] - 10)
subby<-subset(subby,year<= y[12] + 10)
n[12]
a[12]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[12]-10,y[12]+10),type="l",ylab="Polity",xlab="Year",main="Gabon: Age = 65")
abline(v=y[12],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[13])
subby<-subset(subby,year>= y[13] - 10)
subby<-subset(subby,year<= y[13] + 10)
n[13]
a[13]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[13]-10,y[13]+10),type="l",ylab="Polity",xlab="Year",main="Poland: Age = 67")
abline(v=y[13],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
cord.x<-c(dem$year[1],dem$year[1],dem$year[2],dem$year[2])
cord.y<-c(-11,11,11,-11)
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))

subby<-subset(data,ccode==c[14])
subby<-subset(subby,year>= y[14] - 10)
subby<-subset(subby,year<= y[14] + 10)
n[14]
a[14]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[14]-10,y[14]+10),type="l",ylab="Polity",xlab="Year",main="Syria: Age = 70")
abline(v=y[14],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[15])
subby<-subset(subby,year>= y[15] - 10)
subby<-subset(subby,year<= y[15] + 10)
n[15]
a[15]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[15]-10,y[15]+10),type="l",ylab="Polity",xlab="Year",main="Laos: Age = 72")
abline(v=y[15],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[17])
subby<-subset(subby,year>= y[17] - 10)
subby<-subset(subby,year<= y[17] + 10)
n[17]
a[17]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[17]-10,y[17]+10),type="l",ylab="Polity",xlab="Year",main="Bulgaria: Age = 73")
abline(v=y[17],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[19])
subby<-subset(subby,year>= y[19] - 10)
subby<-subset(subby,year<= y[19] + 10)
n[19]
a[19]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[19]-10,y[19]+10),type="l",ylab="Polity",xlab="Year",main="Russia: Age = 74")
abline(v=y[19],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[21])
subby<-subset(subby,year>= y[21] - 10)
subby<-subset(subby,year<= y[20] + 10)
n[20]
a[20]
n[21]
a[21]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[21]-10,y[20]+10),type="l",ylab="Polity",xlab="Year",main="Russia: Age = 76, 74")
abline(v=y[21],lty=1,col="red")	
abline(v=y[20],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[22])
subby<-subset(subby,year>= y[22] - 10)
subby<-subset(subby,year<= y[22] + 10)
n[22]
a[22]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[22]-10,y[22]+10),type="l",ylab="Polity",xlab="Year",main="Liberia: Age = 76")
abline(v=y[22],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[23])
subby<-subset(subby,year>= y[23] - 10)
subby<-subset(subby,year<= y[23] + 10)
n[23]
a[23]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[23]-10,y[23]+10),type="l",ylab="Polity",xlab="Year",main="Albania: Age = 77")
abline(v=y[23],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
cord.x<-c(dem$year[1],dem$year[1],dem$year[5],dem$year[5])
cord.y<-c(-11,11,11,-11)
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))

subby<-subset(data,ccode==c[24])
subby<-subset(subby,year>= y[24] - 10)
subby<-subset(subby,year<= y[24] + 10)
n[24]
a[24]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[24]-10,y[24]+10),type="l",ylab="Polity",xlab="Year",main="Taiwan: Age = 78")
abline(v=y[24],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
cord.x<-c(dem$year[1],dem$year[1],dem$year[3],dem$year[3])
cord.y<-c(-11,11,11,-11)
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))

subby<-subset(data,ccode==c[25])
subby<-subset(subby,year>= y[25] - 10)
subby<-subset(subby,year<= y[25] + 10)
n[25]
a[25]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[25]-10,y[25]+10),type="l",ylab="Polity",xlab="Year",main="Vietnam: Age = 78")
abline(v=y[25],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[26])
subby<-subset(subby,year>= y[26] - 10)
subby<-subset(subby,year<= y[26] + 10)
n[26]
a[26]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[26]-10,y[26]+10),type="l",ylab="Polity",xlab="Year",main="North Korea: Age = 82")
abline(v=y[26],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[27])
subby<-subset(subby,year>= y[27] - 10)
subby<-subset(subby,year<= y[27] + 10)
n[27]
a[27]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[27]-10,y[27]+10),type="l",ylab="Polity",xlab="Year",main="China: Age = 83")
abline(v=y[27],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[28])
subby<-subset(subby,year>= y[28] - 10)
subby<-subset(subby,year<= y[28] + 10)
n[28]
a[28]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[28]-10,y[28]+10),type="l",ylab="Polity",xlab="Year",main="South Africa: Age = 85")
abline(v=y[28],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[29])
subby<-subset(subby,year>= y[29] - 10)
subby<-subset(subby,year<= y[29] + 10)
n[29]
a[29]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[29]-10,y[29]+10),type="l",ylab="Polity",xlab="Year",main="Kenya: Age = 87")
abline(v=y[29],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[30])
subby<-subset(subby,year>= y[30] - 10)
subby<-subset(subby,year<= y[30] + 10)
n[30]
a[30]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[30]-10,y[30]+10),type="l",ylab="Polity",xlab="Year",main="Iran: Age = 87")
abline(v=y[30],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[31])
subby<-subset(subby,year>= y[31] - 10)
subby<-subset(subby,year<= y[31] + 10)
n[31]
a[31]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[31]-10,y[31]+10),type="l",ylab="Polity",xlab="Year",main="Yugoslavia: Age = 88")
abline(v=y[31],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[33])
subby<-subset(subby,year>= y[33] - 10)
subby<-subset(subby,year<= y[33] + 10)
n[33]
a[33]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[33]-10,y[33]+10),type="l",ylab="Polity",xlab="Year",main="Cote d'Ivoire: Age = 88")
abline(v=y[33],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[34])
subby<-subset(subby,year>= y[34] - 10)
subby<-subset(subby,year<= y[18] + 10)
n[34]
a[34]
n[18]
a[18]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[34]-10,y[18]+10),type="l",ylab="Polity",xlab="Year",main="Taiwan: Age = 89, 73")
abline(v=y[34],lty=1,col="red")	
abline(v=y[18],lty=1,col="red")
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[35])
subby<-subset(subby,year>= y[35] - 10)
subby<-subset(subby,year<= y[35] + 10)
n[35]
a[35]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[35]-10,y[35]+10),type="l",ylab="Polity",xlab="Year",main="China: Age = 93")
abline(v=y[35],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

dev.off()






## High Development Countries (Appendix, Figure 5)

c<-high$ccode
y<-high$year
n<-high$ccname
a<-high$agepre

pdf(file="HighDev.pdf",8,8)
par(mfrow=c(4,4))
subby<-subset(data,ccode==c[1])
subby<-subset(subby,year>= y[1] - 10)
subby<-subset(subby,year<= y[1] + 10)
n[1]
a[1]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[1]-10,y[1]+10),type="l",ylab="Polity",xlab="Year",main="Mauritania: Age = 43")
abline(v=y[1],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[2])
subby<-subset(subby,year>= y[2] - 10)
subby<-subset(subby,year<= y[2] + 10)
n[2]
a[2]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[2]-10,y[2]+10),type="l",ylab="Polity",xlab="Year",main="Panama: Age = 52")
abline(v=y[2],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
cord.x<-c(dem$year[1],dem$year[1],dem$year[3],dem$year[3])
cord.y<-c(-11,11,11,-11)
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))


subby<-subset(data,ccode==c[3])
subby<-subset(subby,year>= y[3] - 10)
subby<-subset(subby,year<= y[3] + 10)
n[3]
a[3]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[3]-10,y[3]+10),type="l",ylab="Polity",xlab="Year",main="Kuwait: Age = 64")
abline(v=y[3],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)


subby<-subset(data,ccode==c[4])
subby<-subset(subby,year>= y[4] - 10)
subby<-subset(subby,year<= y[4] + 10)
n[4]
a[4]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[4]-10,y[4]+10),type="l",ylab="Polity",xlab="Year",main="Bahrain: Age = 66")
abline(v=y[4],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[5])
subby<-subset(subby,year>= y[5] - 10)
subby<-subset(subby,year<= y[5] + 10)
n[5]
a[5]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[5]-10,y[5]+10),type="l",ylab="Polity",xlab="Year",main="Poland: Age = 67")
abline(v=y[5],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
cord.x<-c(dem$year[1],dem$year[1],dem$year[2],dem$year[2])
cord.y<-c(-11,11,11,-11)
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))

subby<-subset(data,ccode==c[6])
subby<-subset(subby,year>= y[6] - 10)
subby<-subset(subby,year<= y[6] + 10)
n[6]
a[6]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[6]-10,y[6]+10),type="l",ylab="Polity",xlab="Year",main="Saudi Arabia: Age = 69")
abline(v=y[6],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[7])
subby<-subset(subby,year>= y[7] - 10)
subby<-subset(subby,year<= y[7] + 10)
n[7]
a[7]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[7]-10,y[7]+10),type="l",ylab="Polity",xlab="Year",main="Syria: Age = 70")
abline(v=y[7],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[8])
subby<-subset(subby,year>= y[8] - 10)
subby<-subset(subby,year<= y[8] + 10)
n[8]
a[8]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[8]-10,y[8]+10),type="l",ylab="Polity",xlab="Year",main="Kuwait: Age = 70")
abline(v=y[8],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[9])
subby<-subset(subby,year>= y[9] - 10)
subby<-subset(subby,year<= y[10] + 10)
n[9]
a[9]
a[10]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[9]-10,y[10]+10),type="l",ylab="Polity",xlab="Year",main="Russia: Age = 74,76")
abline(v=y[9],lty=1,col="red")
abline(v=y[10],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
subby<-subset(data,ccode==c[10])
subby<-subset(subby,year<= y[10] + 10)
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[11])
subby<-subset(subby,year>= y[11] - 10)
subby<-subset(subby,year<= y[11] + 10)
n[11]
a[11]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[11]-10,y[11]+10),type="l",ylab="Polity",xlab="Year",main="Liberia: Age = 76")
abline(v=y[11],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[12])
subby<-subset(subby,year>= y[12] - 10)
subby<-subset(subby,year<= y[12] + 10)
n[12]
a[12]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[12]-10,y[12]+10),type="l",ylab="Polity",xlab="Year",main="Taiwan: Age = 78")
abline(v=y[12],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
cord.x<-c(dem$year[1],dem$year[1],dem$year[3],dem$year[3])
cord.y<-c(-11,11,11,-11)
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))

subby<-subset(data,ccode==c[13])
subby<-subset(subby,year>= y[13] - 10)
subby<-subset(subby,year<= y[13] + 10)
n[13]
a[13]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[13]-10,y[13]+10),type="l",ylab="Polity",xlab="Year",main="Portugal: Age = 79")
abline(v=y[13],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
cord.x<-c(dem$year[1],dem$year[1],dem$year[3],dem$year[3])
cord.y<-c(-11,11,11,-11)
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))

subby<-subset(data,ccode==c[14])
subby<-subset(subby,year>= y[14] - 10)
subby<-subset(subby,year<= y[14] + 10)
n[14]
a[14]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[14]-10,y[14]+10),type="l",ylab="Polity",xlab="Year",main="Spain: Age = 83")
abline(v=y[14],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
cord.x<-c(dem$year[1],dem$year[1],dem$year[9],dem$year[9])
cord.y<-c(-11,11,11,-11)
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))

dev.off()


## Low Development Countries (Appendix, Figure 6)


c<-low$ccode
y<-low$year
n<-low$ccname
a<-low$agepre

pdf(file="LowDev.pdf",8,8)
par(mfrow=c(4,4))
subby<-subset(data,ccode==c[1])
subby<-subset(subby,year>= y[1] - 10)
subby<-subset(subby,year<= y[1] + 10)
n[1]
a[1]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[1]-10,y[1]+10),type="l",ylab="Polity",xlab="Year",main="Nepal: Age = 49")
abline(v=y[1],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[2])
subby<-subset(subby,year>= y[2] - 10)
subby<-subset(subby,year<= y[2] + 10)
n[2]
a[2]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[2]-10,y[2]+10),type="l",ylab="Polity",xlab="Year",main="Nepal: Age = 52")
abline(v=y[2],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[3])
subby<-subset(subby,year>= y[3] - 10)
subby<-subset(subby,year<= y[3] + 10)
n[3]
a[3]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[3]-10,y[3]+10),type="l",ylab="Polity",xlab="Year",main="Nigeria: Age = 55")
abline(v=y[3],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
cord.x<-c(dem$year[1],dem$year[1],dem$year[10],dem$year[10])
cord.y<-c(-11,11,11,-11)
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))

subby<-subset(data,ccode==c[4])
subby<-subset(subby,year>= y[4] - 10)
subby<-subset(subby,year<= y[4] + 10)
n[4]
a[4]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[4]-10,y[4]+10),type="l",ylab="Polity",xlab="Year",main="Mongolia: Age = 57")
abline(v=y[4],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[5])
subby<-subset(subby,year>= y[5] - 10)
subby<-subset(subby,year<= y[5] + 10)
n[5]
a[5]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[5]-10,y[5]+10),type="l",ylab="Polity",xlab="Year",main="Guinea: Age = 62")
abline(v=y[5],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[6])
subby<-subset(subby,year>= y[6] - 10)
subby<-subset(subby,year<= y[6] + 10)
n[6]
a[6]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[6]-10,y[6]+10),type="l",ylab="Polity",xlab="Year",main="Comoros: Age = 62")
abline(v=y[6],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)
cord.x<-c(dem$year[1],dem$year[1],dem$year[5],dem$year[5])
cord.y<-c(-11,11,11,-11)
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))
cord.x<-c(dem$year[6],dem$year[6],dem$year[10],dem$year[10])
cord.y<-c(-11,11,11,-11)
polygon(cord.x,cord.y,col=rgb(0,0,0,alpha=0.3))

subby<-subset(data,ccode==c[7])
subby<-subset(subby,year>= y[7] - 10)
subby<-subset(subby,year<= y[7] + 10)
n[7]
a[7]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[7]-10,y[7]+10),type="l",ylab="Polity",xlab="Year",main="Haiti: Age = 64")
abline(v=y[7],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[8])
subby<-subset(subby,year>= y[8] - 10)
subby<-subset(subby,year<= y[8] + 10)
n[8]
a[8]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[8]-10,y[8]+10),type="l",ylab="Polity",xlab="Year",main="Gabon: Age = 65")
abline(v=y[8],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)


subby<-subset(data,ccode==c[9])
subby<-subset(subby,year>= y[9] - 10)
subby<-subset(subby,year<= y[9] + 10)
n[9]
a[9]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[9]-10,y[9]+10),type="l",ylab="Polity",xlab="Year",main="Laos: Age = 72")
abline(v=y[9],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[10])
subby<-subset(subby,year>= y[10] - 10)
subby<-subset(subby,year<= y[10] + 10)
n[10]
a[10]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[10]-10,y[10]+10),type="l",ylab="Polity",xlab="Year",main="Vietnam: Age = 78")
abline(v=y[10],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[11])
subby<-subset(subby,year>= y[11] - 10)
subby<-subset(subby,year<= y[11] + 10)
n[11]
a[11]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[11]-10,y[11]+10),type="l",ylab="Polity",xlab="Year",main="Swaziland: Age = 83")
abline(v=y[11],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[12])
subby<-subset(subby,year>= y[12] - 10)
subby<-subset(subby,year<= y[12] + 10)
n[12]
a[12]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[12]-10,y[12]+10),type="l",ylab="Polity",xlab="Year",main="China: Age = 83")
abline(v=y[12],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[13])
subby<-subset(subby,year>= y[13] - 10)
subby<-subset(subby,year<= y[13] + 10)
n[13]
a[13]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[13]-10,y[13]+10),type="l",ylab="Polity",xlab="Year",main="South Africa: Age = 85")
abline(v=y[13],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

subby<-subset(data,ccode==c[14])
subby<-subset(subby,year>= y[14] - 10)
subby<-subset(subby,year<= y[14] + 10)
n[14]
a[14]
plot(subby$year,subby$polity2, ylim=c(-10,10),xlim=c(y[14]-10,y[14]+10),type="l",ylab="Polity",xlab="Year",main="Kenya: Age = 87")
abline(v=y[14],lty=1,col="red")	
dem<-subset(subby,democracy==1)
nrow(dem)

dev.off()



#Ultimate fate of dictators, by regime type (Figure 5 in Main Article)

nondem<-subset(data,polity_lag1<0)
nondem<-subset(nondem, is.na(death)==FALSE )
subby<-nondem
subby$keep<-subby$otherexit
subby$keep<-ifelse(subby$death==1,1,subby$keep)
subby<-subset(subby,keep==1)


#subby<-merge(subby,leaders,all.x=TRUE,all.y=FALSE)

death<-subset(nondem,death==1)
age<-quantile(death$agepre,c(0.25,0.5,0.75),na.rm=TRUE)

age0<-subset(subby,agepre<=54)

age25<-subset(subby,agepre>54)
age25<-subset(age25,agepre<=64)

age50<-subset(subby,subby$agepre>64)
age50<-subset(age50,age50$agepre<=73)

age75<-subset(subby,agepre>73)


#non personalist

age0x<-subset(age0,personal==0)
age25x<-subset(age25,personal==0)
age50x<-subset(age50,personal==0)
age75x<-subset(age75,personal==0)

death0<-nrow(subset(age0x,death==1))
s1<-subset(age0x,democracy==1)
s2<-subset(age0x,democracy==0)
irreg0<-nrow(subset(s2,otherirreg==1))
reg0<-subset(s2,otherexit==1)
reg0<-nrow(subset(reg0,otherirreg==0))
dem0<-nrow(subset(s1,death==0))

death25<-nrow(subset(age25x,death==1))
s1<-subset(age25x,democracy==1)
s2<-subset(age25x,democracy==0)
irreg25<-nrow(subset(s2,otherirreg==1))
reg25<-subset(s2,otherexit==1)
reg25<-nrow(subset(reg25,otherirreg==0))
dem25<-nrow(subset(s1,death==0))

death50<-nrow(subset(age50x,death==1))
s1<-subset(age50x,democracy==1)
s2<-subset(age50x,democracy==0)
irreg50<-nrow(subset(s2,otherirreg==1))
reg50<-subset(s2,otherexit==1)
reg50<-nrow(subset(reg50,otherirreg==0))
dem50<-nrow(subset(s1,death==0))

death75<-nrow(subset(age75x,death==1))
s1<-subset(age75x,democracy==1)
s2<-subset(age75x,democracy==0)
irreg75<-nrow(subset(s2,otherirreg==1))
reg75<-subset(s2,otherexit==1)
reg75<-nrow(subset(reg75,otherirreg==0))
dem75<-nrow(subset(s1,death==0))

death<-c(death0,death25,death50,death75)
irreg<-c(irreg0,irreg25,irreg50,irreg75)
reg<-c(reg0,reg25,reg50,reg75)
dem<-c(dem0,dem25,dem50,dem75)

noper_exit<-data.frame(death,irreg,reg,dem)


# personal
age0x<-subset(age0,personal==1)
age25x<-subset(age25,personal==1)
age50x<-subset(age50,personal==1)
age75x<-subset(age75,personal==1)

death0<-nrow(subset(age0x,death==1))
s1<-subset(age0x,democracy==1)
s2<-subset(age0x,democracy==0)
irreg0<-nrow(subset(s2,otherirreg==1))
reg0<-subset(s2,otherexit==1)
reg0<-nrow(subset(reg0,otherirreg==0))
dem0<-nrow(subset(s1,death==0))

death25<-nrow(subset(age25x,death==1))
s1<-subset(age25x,democracy==1)
s2<-subset(age25x,democracy==0)
irreg25<-nrow(subset(s2,otherirreg==1))
reg25<-subset(s2,otherexit==1)
reg25<-nrow(subset(reg25,otherirreg==0))
dem25<-nrow(subset(s1,death==0))

death50<-nrow(subset(age50x,death==1))
s1<-subset(age50x,democracy==1)
s2<-subset(age50x,democracy==0)
irreg50<-nrow(subset(s2,otherirreg==1))
reg50<-subset(s2,otherexit==1)
reg50<-nrow(subset(reg50,otherirreg==0))
dem50<-nrow(subset(s1,death==0))

death75<-nrow(subset(age75x,death==1))
s1<-subset(age75x,democracy==1)
s2<-subset(age75x,democracy==0)
irreg75<-nrow(subset(s2,otherirreg==1))
reg75<-subset(s2,otherexit==1)
reg75<-nrow(subset(reg75,otherirreg==0))
dem75<-nrow(subset(s1,death==0))

death<-c(death0,death25,death50,death75)
irreg<-c(irreg0,irreg25,irreg50,irreg75)
reg<-c(reg0,reg25,reg50,reg75)
dem<-c(dem0,dem25,dem50,dem75)

per_exit<-data.frame(death,irreg,reg,dem)


y<-per_exit

x0<-sum(y[1,])+sum(y[2,])+sum(y[3,])+sum(y[4,])
x25<-sum(y[2,])+sum(y[3,])+sum(y[4,])
x50<-sum(y[3,])+sum(y[4,])
x75<-sum(y[4,])

a0<-c(sum(y[,1])/x0,sum(y[,2])/x0,sum(y[,3])/x0,sum(y[,4])/x0)
y<-y[-1,]
a25<-c(sum(y[,1])/x25,sum(y[,2])/x25,sum(y[,3])/x25,sum(y[,4])/x25)
y<-y[-1,]
a50<-c(sum(y[,1])/x50,sum(y[,2])/x50,sum(y[,3])/x50,sum(y[,4])/x50)
y<-y[-1,]
a75<-c(sum(y[,1])/x75,sum(y[,2])/x75,sum(y[,3])/x75,sum(y[,4])/x75)

per_perc<-rbind(a0,a25,a50,a75)
colnames(per_perc)<-c("death","irreg","reg","dem")

y<-noper_exit

x0<-sum(y[1,])+sum(y[2,])+sum(y[3,])+sum(y[4,])
x25<-sum(y[2,])+sum(y[3,])+sum(y[4,])
x50<-sum(y[3,])+sum(y[4,])
x75<-sum(y[4,])

a0<-c(sum(y[,1])/x0,sum(y[,2])/x0,sum(y[,3])/x0,sum(y[,4])/x0)
y<-y[-1,]
a25<-c(sum(y[,1])/x25,sum(y[,2])/x25,sum(y[,3])/x25,sum(y[,4])/x25)
y<-y[-1,]
a50<-c(sum(y[,1])/x50,sum(y[,2])/x50,sum(y[,3])/x50,sum(y[,4])/x50)
y<-y[-1,]
a75<-c(sum(y[,1])/x75,sum(y[,2])/x75,sum(y[,3])/x75,sum(y[,4])/x75)

noper_perc<-rbind(a0,a25,a50,a75)
colnames(noper_perc)<-c("death","irreg","reg","dem")

personalist<-c(per_perc[3,1],per_perc[3,2],per_perc[3,3],per_perc[3,4])
nonpersonalist<-c(noper_perc[3,1],noper_perc[3,2],noper_perc[3,3],noper_perc[3,4])

plotty<-data.frame(personalist, nonpersonalist)
plotty<-as.matrix(plotty)

pdf(file="fig_exits.pdf")
barplot(height=plotty,col=c("black","darkgrey","lightgrey","white"), beside=T, legend=c("death","irreg","reg","dem"), ylim=c(0,0.6))
dev.off()


#Ultimate fate of dictators, by development level (Figure 9 in Main Article)

nondem<-subset(data,polity_lag1<0)

q1<-quantile(nondem$gdppcln_lag1,.75,na.rm=TRUE) 
q2<-quantile(nondem$gdppcln_lag1,.25,na.rm=TRUE) 


#high development

age0x<-subset(age0,gdppcln_lag1>q1)
age25x<-subset(age25,gdppcln_lag1>q1)
age50x<-subset(age50,gdppcln_lag1>q1)
age75x<-subset(age75,gdppcln_lag1>q1)

death0<-nrow(subset(age0x,death==1))
s1<-subset(age0x,democracy==1)
s2<-subset(age0x,democracy==0)
irreg0<-nrow(subset(s2,otherirreg==1))
reg0<-subset(s2,otherexit==1)
reg0<-nrow(subset(reg0,otherirreg==0))
dem0<-nrow(subset(s1,death==0))

death25<-nrow(subset(age25x,death==1))
s1<-subset(age25x,democracy==1)
s2<-subset(age25x,democracy==0)
irreg25<-nrow(subset(s2,otherirreg==1))
reg25<-subset(s2,otherexit==1)
reg25<-nrow(subset(reg25,otherirreg==0))
dem25<-nrow(subset(s1,death==0))

death50<-nrow(subset(age50x,death==1))
s1<-subset(age50x,democracy==1)
s2<-subset(age50x,democracy==0)
irreg50<-nrow(subset(s2,otherirreg==1))
reg50<-subset(s2,otherexit==1)
reg50<-nrow(subset(reg50,otherirreg==0))
dem50<-nrow(subset(s1,death==0))

death75<-nrow(subset(age75x,death==1))
s1<-subset(age75x,democracy==1)
s2<-subset(age75x,democracy==0)
irreg75<-nrow(subset(s2,otherirreg==1))
reg75<-subset(s2,otherexit==1)
reg75<-nrow(subset(reg75,otherirreg==0))
dem75<-nrow(subset(s1,death==0))

death<-c(death0,death25,death50,death75)
irreg<-c(irreg0,irreg25,irreg50,irreg75)
reg<-c(reg0,reg25,reg50,reg75)
dem<-c(dem0,dem25,dem50,dem75)

high_exit<-data.frame(death,irreg,reg,dem)


#low development

age0x<-subset(age0,gdppcln_lag1<q1)
age25x<-subset(age25,gdppcln_lag1<q1)
age50x<-subset(age50,gdppcln_lag1<q1)
age75x<-subset(age75,gdppcln_lag1<q1)

death0<-nrow(subset(age0x,death==1))
s1<-subset(age0x,democracy==1)
s2<-subset(age0x,democracy==0)
irreg0<-nrow(subset(s2,otherirreg==1))
reg0<-subset(s2,otherexit==1)
reg0<-nrow(subset(reg0,otherirreg==0))
dem0<-nrow(subset(s1,death==0))

death25<-nrow(subset(age25x,death==1))
s1<-subset(age25x,democracy==1)
s2<-subset(age25x,democracy==0)
irreg25<-nrow(subset(s2,otherirreg==1))
reg25<-subset(s2,otherexit==1)
reg25<-nrow(subset(reg25,otherirreg==0))
dem25<-nrow(subset(s1,death==0))

death50<-nrow(subset(age50x,death==1))
s1<-subset(age50x,democracy==1)
s2<-subset(age50x,democracy==0)
irreg50<-nrow(subset(s2,otherirreg==1))
reg50<-subset(s2,otherexit==1)
reg50<-nrow(subset(reg50,otherirreg==0))
dem50<-nrow(subset(s1,death==0))

death75<-nrow(subset(age75x,death==1))
s1<-subset(age75x,democracy==1)
s2<-subset(age75x,democracy==0)
irreg75<-nrow(subset(s2,otherirreg==1))
reg75<-subset(s2,otherexit==1)
reg75<-nrow(subset(reg75,otherirreg==0))
dem75<-nrow(subset(s1,death==0))

death<-c(death0,death25,death50,death75)
irreg<-c(irreg0,irreg25,irreg50,irreg75)
reg<-c(reg0,reg25,reg50,reg75)
dem<-c(dem0,dem25,dem50,dem75)

low_exit<-data.frame(death,irreg,reg,dem)


y<-high_exit

x0<-sum(y[1,])+sum(y[2,])+sum(y[3,])+sum(y[4,])
x25<-sum(y[2,])+sum(y[3,])+sum(y[4,])
x50<-sum(y[3,])+sum(y[4,])
x75<-sum(y[4,])

a0<-c(sum(y[,1])/x0,sum(y[,2])/x0,sum(y[,3])/x0,sum(y[,4])/x0)
y<-y[-1,]
a25<-c(sum(y[,1])/x25,sum(y[,2])/x25,sum(y[,3])/x25,sum(y[,4])/x25)
y<-y[-1,]
a50<-c(sum(y[,1])/x50,sum(y[,2])/x50,sum(y[,3])/x50,sum(y[,4])/x50)
y<-y[-1,]
a75<-c(sum(y[,1])/x75,sum(y[,2])/x75,sum(y[,3])/x75,sum(y[,4])/x75)

high_perc<-rbind(a0,a25,a50,a75)
colnames(high_perc)<-c("death","irreg","reg","dem")


y<-low_exit

x0<-sum(y[1,])+sum(y[2,])+sum(y[3,])+sum(y[4,])
x25<-sum(y[2,])+sum(y[3,])+sum(y[4,])
x50<-sum(y[3,])+sum(y[4,])
x75<-sum(y[4,])

a0<-c(sum(y[,1])/x0,sum(y[,2])/x0,sum(y[,3])/x0,sum(y[,4])/x0)
y<-y[-1,]
a25<-c(sum(y[,1])/x25,sum(y[,2])/x25,sum(y[,3])/x25,sum(y[,4])/x25)
y<-y[-1,]
a50<-c(sum(y[,1])/x50,sum(y[,2])/x50,sum(y[,3])/x50,sum(y[,4])/x50)
y<-y[-1,]
a75<-c(sum(y[,1])/x75,sum(y[,2])/x75,sum(y[,3])/x75,sum(y[,4])/x75)

low_perc<-rbind(a0,a25,a50,a75)
colnames(low_perc)<-c("death","irreg","reg","dem")


low<-c(low_perc[3,1],low_perc[3,2],low_perc[3,3],low_perc[3,4])
high<-c(high_perc[3,1],high_perc[3,2],high_perc[3,3],high_perc[3,4])

plotty<-data.frame(high,low)
plotty<-as.matrix(plotty)

pdf(file="fig_exitsEcon.pdf")
barplot(height=plotty,col=c("black","darkgrey","lightgrey","white"), beside=T, legend=c("death","irreg","reg","dem"), ylim=c(0,0.6),names.arg=c("high development", "low development"))
dev.off()

### Pre- and Post-Death Polity Scores (Figure 1 in the Main Article)

death<-subset(data,death==1)

pdf(file="fig_prepostpolity1.pdf")
plot(death$polity_lag1,death$polity_for1,xlim=c(-10,10),ylim=c(-10,10),ylab="Polity, 1 year after",xlab="Polity, year before")
abline(a=0,b=1)
abline(a=5,b=1,lty=5)
abline(a=-5,b=1,lty=5)
s<-subset(death, abs(polity_for1-polity_lag1)>5)

text(8,-8.5,"Greece 1936",cex=0.75) 
text(4,-5.5,"Ethiopia 1930",cex=0.75) 
text(4.2,-2.5,"Comoros 1998",cex=0.75) 
text(2,-9.5,"Paraguay 1940",cex=0.75) 

text(-7,1.5,"Spain 1975",cex=0.75) 
text(-5,8.5,"Croatia 1999",cex=0.75) 
text(-9,0.5,"Greece 1941",cex=0.75) 
text(-9,-1.5,"Bulgaria 1943",cex=0.75) 
text(-6,4.5,"Nigeria 1998",cex=0.75) 
text(-6.1,0.5,"Yemen 1962",cex=0.75) 
dev.off()

pdf(file="fig_prepostpolity5.pdf")
plot(death$polity_lag1,death$polity_for5,xlim=c(-10,10),ylim=c(-10,10),ylab="Ave. Polity, 5 years after",xlab="Polity, year before")
abline(a=0,b=1)
abline(a=5,b=1,lty=5)
abline(a=-5,b=1,lty=5)
s<-subset(death, abs(polity_for5-polity_lag1)>5)
text(-7,7.1,"Spain 1975",cex=0.75) 
text(-5,8.5,"Croatia 1999",cex=0.75) 
text(-8,5.3,"Greece 1941",cex=0.75) 
text(-6,4.5,"Nigeria 1998",cex=0.75) 

text(8,-1.3,"Colombia 1884",cex=0.75) 
text(2,-9.5,"Paraguay 1940",cex=0.75) 
text(6,-5.5,"Argentina 1974",cex=0.75) 
text(8,-6.9,"Greece 1936",cex=0.75) 
text(4,-4.5,"Ethiopia 1930",cex=0.75) 
dev.off()


pdf(file="fig_prepostpolity10.pdf")
plot(death$polity_lag1,death$polity_for10,xlim=c(-10,10),ylim=c(-10,10),ylab="Ave. Polity, 10 years after",xlab="Polity, year before")
abline(a=0,b=1)
abline(a=5,b=1,lty=5)
abline(a=-5,b=1,lty=5)
s<-subset(death, abs(polity_for10-polity_lag1)>5)

text(2,-5.1,"Guatemala 1926",cex=0.6) 
text(8,0.8,"Colombia 1882",cex=0.6) 
text(8,-2.4,"Colombia 1884",cex=0.6) 
text(2,-8,"Paraguay 1940",cex=0.6) 
text(6,-4.4,"Argentina 1974",cex=0.6) 
text(8.6,2.1,"Uruguay 1965",cex=0.6) 
text(8,-1.1,"Uruguay 1967",cex=0.6) 
text(8,-0.3,"Greece 1936",cex=0.6) 
text(5.6,-1.9,"Sierra Leone 1964",cex=0.6) 

text(-7.6,-1.4,"Guyana 1985",cex=0.6) 
text(-7,8.7,"Spain 1975",cex=0.6) 
text(-9.5,-2.3,"Portugal 1968",cex=0.6) 
text(-9.5,-0.9,"Albania 1985",cex=0.6) 
text(-5,9,"Croatia 1999",cex=0.6) 
text(-8,5.7,"Greece 1941",cex=0.6) 
text(-7,-0.5,"Russia 1985",cex=0.6) 
text(-7,0.5,"Niger 1987",cex=0.6) 
text(-6,4.5,"Nigeria 1998",cex=0.6) 
text(-1,5.6,"Taiwan",cex=0.6) 
dev.off()


##Treisman replications

library(MASS)
library(foreign)

library(survival)
library(stargazer)
library(lmtest)

library(sandwich)

death<-read.csv("data_TreismanDemocworki.csv")

death$ncfull<-ifelse((
death$country=="Nepal" & death$year==1877|
death$country=="Mongolia" & death$year==1952|
death$country=="Nepal" & death$year==1901| 
death$country=="Nepal" & death$year==1948|
death$country=="Nepal" & death$year==1932|
death$country=="Nepal" & death$year==1929|
death$country=="Guinea" & death$year==1984|
death$country=="Nepal" & death$year==1955| 
death$country=="China" & death$year==1908|
death$country=="China" & death$year==1916|
death$country=="Niger" & death$year==1987|
death$country=="Comoros" & death$year==1998|
death$country=="Venezuela" & death$year==1878|
death$country=="Nepal" & death$year==1972|
death$country=="Morocco" & death$year==1894|
death$country=="Brazil" & death$year==1909|
death$country=="Peru" & death$year==1904|
death$country=="Thailand" & death$year==1925|
death$country=="Thailand" & death$year==1910|
death$country=="China" & death$year==1976|
death$country=="Philippines" & death$year==1948|
death$country=="Egypt" & death$year==1936|
death$country=="Haiti" & death$year==1971|
death$country=="Mozambique" & death$year==1986|
death$country=="Vietnam" & death$year==1986|
death$country=="Laos" & death$year==1992|
death$country=="Iran" & death$year==1907|
death$country=="Angola" & death$year==1979|
death$country=="Mauritania" & death$year==1979|
death$country=="Yemen North" & death$year==1962|
death$country=="Kenya" & death$year==1978|
death$country=="Portugal" & death$year==1889|
death$country=="Thailand" & death$year==1963|
death$country=="Nigeria" & death$year==1998|
death$country=="Russia" & death$year==1894|
death$country=="Nicaragua" & death$year==1923|
death$country=="Iraq" & death$year==1933|
death$country=="Egypt" & death$year==1970|
death$country=="Ivory Coast" & death$year==1993|
death$country=="Iraq" & death$year==1939|
death$country=="Romania" & death$year==1927|
death$country=="Morocco" & death$year==1961|
death$country=="Philippines" & death$year==1957|
death$country=="USSR" & death$year==1923|
death$country=="Bulgaria" & death$year==1949|
death$country=="Guatemala" & death$year==1926|
death$country=="Bulgaria" & death$year==1943|
death$country=="Liberia" & death$year==1971|
death$country=="Bulgaria" & death$year==1950|
death$country=="Poland" & death$year==1935|
death$country=="Turkey" & death$year==1938|
death$country=="Italy" & death$year==1887|
death$country=="South Africa" & death$year==1919|
death$country=="Romania" & death$year==1914|
death$country=="Japan" & death$year==1923|
death$country=="Japan" & death$year==1926|
death$country=="Paraguay" & death$year==1940|
death$country=="Bolivia" & death$year==1969|
death$country=="Greece" & death$year==1941|
death$country=="Romania" & death$year==1965|
death$country=="Germany" & death$year==1888|
death$country=="Greece" & death$year==1955|
death$country=="Albania" & death$year==1985|
death$country=="Saudi Arabia" & death$year==1953|
death$country=="Swaziland" & death$year==1982|
death$country=="Korea North" & death$year==1994|
death$country=="Malaysia" & death$year==1976|
death$country=="Austria" & death$year==1916|
death$country=="Morocco" & death$year==1999|
death$country=="Chile" & death$year==1910|
death$country=="Nicaragua" & death$year==1966|
death$country=="Algeria" & death$year==1978|
death$country=="Poland" & death$year==1956|
death$country=="China" & death$year==1997|
death$country=="USSR" & death$year==1953|
death$country=="South Africa" & death$year==1958|
death$country=="Venezuela" & death$year==1935|
death$country=="Chile" & death$year==1941|
death$country=="Iraq" & death$year==1966|
death$country=="Iran" & death$year==1989|
death$country=="Taiwan" & death$year==1975|
death$country=="Argentina" & death$year==1906|
death$country=="Czechoslovakia" & death$year==1953|
death$country=="Jordan" & death$year==1999|
death$country=="Czechoslovakia" & death$year==1957|
death$country=="Taiwan" & death$year==1978|
death$country=="Portugal" & death$year==1968|
death$country=="Bahrain" & death$year==1999|
death$country=="Gabon" & death$year==1967|
death$country=="Panama" & death$year==1981|
death$country=="Yugoslavia" & death$year==1980|
death$country=="Poland" & death$year==1980|
death$country=="Croatia" & death$year==1999|
death$country=="USSR" & death$year==1982|
death$country=="USSR" & death$year==1984|
death$country=="USSR" & death$year==1985|
death$country=="Syria" & death$year==2000|
death$country=="Spain" & death$year==1975|
death$country=="Taiwan" & death$year==1988|
death$country=="Saudi Arabia" & death$year==1982|
death$country=="Kuwait" & death$year==1978|
death$country=="Kuwait" & death$year==1965|
death$country=="Haiti" & death$year==1896|
death$country=="Haiti" & death$year==1913|
death$country=="Nicaragua" & death$year==1889|
death$country=="Panama" & death$year==1910|
death$country=="Panama" & death$year==1918|
death$country=="Panama" & death$year==1939|
death$country=="Guyana" & death$year==1985|
death$country=="Ecuador" & death$year==1911|
death$country=="Ecuador" & death$year==1939|
death$country=="Peru" & death$year==1894|
death$country=="Paraguay" & death$year==1880|
death$country=="Paraguay" & death$year==1919|
death$country=="Liberia" & death$year==1896|
death$country=="Ethiopia" & death$year==1911|
death$country=="Ethiopia" & death$year==1930|
death$country=="Orange Free State" & death$year==1888|
death$country=="Oman" & death$year==1888|
death$country=="Oman" & death$year==1913|
death$country=="Afghanistan" & death$year==1880|
death$country=="Afghanistan" & death$year==1901|
death$country=="Bhutan" & death$year==1952|
death$country=="Bhutan" & death$year==1972|
death$country=="Pakistan" & death$year==1948),1,0)

#leaderturn=102
for(i in c(2:nrow(death))){
	k<-ncol(death)
	j<-i-1
	death[i,k]<-ifelse(death[i,3]!=death[j,3] | is.na(death[j,102])==TRUE,NA,death[i,k])
}
death[1,k]<-NA

# pol2norm is 37
death$f10pol2norm <- rep(NA,nrow(death))
for(i in c(1:(nrow(death)-10))){
	k<-ncol(death)
	j<-i+10
	death[i,k]<-ifelse(death[i,3]==death[j,3],death[j,37],NA)
}

death$df10pol2norm<-death$f10pol2norm-death$pol2norm

death$loggdppc <-log10(death$gdppc)
death$lloggdppc<-rep(NA,nrow(death))
for(i in c(2:nrow(death))){
	j<-i-1
	k<-ncol(death)
	l<-k-1
	death[i,k]<-ifelse(death[i,3]==death[j,3],death[j,l],NA)
}


#ncfull is 553
death$nonatc10<-rep(NA, nrow(death))
for(i in c(1:(nrow(death)-10))){
	k<-ncol(death)
	j<-i+10
death[i,k]<-ifelse(sum(death$ncfull[i:j])>0,0,1)
}


#ncfull is 553
death$nolt10<-rep(NA, nrow(death))
for(i in c(1:(nrow(death)-10))){
	k<-ncol(death)
	j<-i+10
death[i,k]<-ifelse(sum(death$leaderturn[i:j])>0,0,1)
}


death$nonatc10_alt<-rep(NA, nrow(death))
for(i in c(1:(nrow(death)-10))){
	k<-ncol(death)
	j<-i+10
	l<-i+1
death[i,k]<-ifelse(sum(death$ncfull[l:j])>0,0,1)
}

death$nolt10_alt<-rep(NA, nrow(death))
for(i in c(1:(nrow(death)-10))){
	k<-ncol(death)
	j<-i+10
	l<-i+1
death[i,k]<-ifelse(sum(death$leaderturn[l:j])>0,0,1)
}


#TRIESMAN MODELS + CORRECT PAIRINGS (Table 33 in online appendix)
all1<-subset(death,(ncfull==1&polity2<6&year>1874&year<2005))
all0<-subset(death,(ncfull==0&polity2<6&year>1874&year<2005))

m1_death<-lm(df10pol2norm~lloggdppc, data=all1)
G<-length(unique(all1$country))
N<-length(all1$country)
dfa<-(G/(G-1))*(N-1)/m1_death$df.residual
m1_death_cluster<-dfa*vcovHC(m1_death,type="HC0",cluster="group",adjust=T)
m1_death_c<-coeftest(m1_death,vcov=m1_death_cluster)

m1_nodeath<-lm(df10pol2norm~lloggdppc, data=all0)
G<-length(unique(all0$country))
N<-length(all0$country)
dfa<-(G/(G-1))*(N-1)/m1_nodeath$df.residual
m1_nodeath_cluster<-dfa*vcovHC(m1_nodeath,type="HC0",cluster="group",adjust=T)
m1_nodeath_c<-coeftest(m1_nodeath,vcov=m1_nodeath_cluster)


nodeath1<-subset(death,(nonatc10_alt==1&ncfull==1&polity2<6&year>1874&year<2005))
nodeath0<-subset(death,(nonatc10==1&ncfull==0&polity2<6&year>1874&year<2005))

m2_death<-lm(df10pol2norm~lloggdppc, data=nodeath1)
G<-length(unique(nodeath1$country))
N<-length(nodeath1$country)
dfa<-(G/(G-1))*(N-1)/m2_death$df.residual
m2_death_cluster<-dfa*vcovHC(m2_death,type="HC0",cluster="group",adjust=T)
m2_death_c<-coeftest(m2_death,vcov=m2_death_cluster)

m2_nodeath<-lm(df10pol2norm~lloggdppc, data=nodeath0)
G<-length(unique(nodeath0$country))
N<-length(nodeath0$country)
dfa<-(G/(G-1))*(N-1)/m2_nodeath$df.residual
m2_nodeath_cluster<-dfa*vcovHC(m2_nodeath,type="HC0",cluster="group",adjust=T)
m2_nodeath_c<-coeftest(m2_nodeath,vcov=m2_nodeath_cluster)


noturn1<-subset(death,(nolt10_alt==1&ncfull==1&polity2<6&year>1874&year<2005))
noturn0<-subset(death,(nolt10==1&ncfull==0&polity2<6&year>1874&year<2005))

m3_death<-lm(df10pol2norm~lloggdppc, data=noturn1)
G<-length(unique(noturn1$country))
N<-length(noturn1$country)
dfa<-(G/(G-1))*(N-1)/m3_death$df.residual
m3_death_cluster<-dfa*vcovHC(m3_death,type="HC0",cluster="group",adjust=T)
m3_death_c<-coeftest(m3_death,vcov=m3_death_cluster)

m3_nodeath<-lm(df10pol2norm~lloggdppc, data=noturn0)
G<-length(unique(noturn0$country))
N<-length(noturn0$country)
dfa<-(G/(G-1))*(N-1)/m3_nodeath$df.residual
m3_nodeath_cluster<-dfa*vcovHC(m3_nodeath,type="HC0",cluster="group",adjust=T)
m3_nodeath_c<-coeftest(m3_nodeath,vcov=m3_nodeath_cluster)






#ALTERNATIVE WAY OF MODELING  (Table 34 in online appendix)


all<-subset(death,polity2<6&year>1874&year<2005)
nodeath<-subset(death,(nonatc10_alt==1&polity2<6&year>1874&year<2005))
noturn<-subset(death,(nolt10_alt==1&polity2<6&year>1874&year<2005))


m1<-lm(df10pol2norm~lloggdppc*ncfull, data=all)
G<-length(unique(all$country))
N<-length(all$country)
dfa<-(G/(G-1))*(N-1)/m1$df.residual
m1_cluster<-dfa*vcovHC(m1,type="HC0",cluster="group",adjust=T)
m1_c<-coeftest(m1,vcov=m1_cluster)


m2<-lm(df10pol2norm~lloggdppc*ncfull, data=nodeath)
G<-length(unique(nodeath$country))
N<-length(nodeath$country)
dfa<-(G/(G-1))*(N-1)/m2$df.residual
m2_cluster<-dfa*vcovHC(m2,type="HC0",cluster="group",adjust=T)
m2_c<-coeftest(m2,vcov=m2_cluster)

m3<-lm(df10pol2norm~lloggdppc*ncfull, data=noturn)
G<-length(unique(noturn$country))
N<-length(noturn$country)
dfa<-(G/(G-1))*(N-1)/m3$df.residual
m3_cluster<-dfa*vcovHC(m3,type="HC0",cluster="group",adjust=T)
m3_c<-coeftest(m3,vcov=m3_cluster)

m4<-lm(df10pol2norm~lloggdppc*ncfull+factor(country), data=all)
G<-length(unique(all$country))
N<-length(all$country)
dfa<-(G/(G-1))*(N-1)/m4$df.residual
m4_cluster<-dfa*vcovHC(m4,type="HC0",cluster="group",adjust=T)
m4_c<-coeftest(m4,vcov=m4_cluster)


m5<-lm(df10pol2norm~lloggdppc*ncfull +factor(country), data=nodeath)
G<-length(unique(nodeath$country))
N<-length(nodeath$country)
dfa<-(G/(G-1))*(N-1)/m5$df.residual
m5_cluster<-dfa*vcovHC(m5,type="HC0",cluster="group",adjust=T)
m5_c<-coeftest(m5,vcov=m5_cluster)

m6<-lm(df10pol2norm~lloggdppc*ncfull+factor(country), data=noturn)
G<-length(unique(noturn$country))
N<-length(noturn$country)
dfa<-(G/(G-1))*(N-1)/m6$df.residual
m6_cluster<-dfa*vcovHC(m6,type="HC0",cluster="group",adjust=T)
m6_c<-coeftest(m6,vcov=m6_cluster)


